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ABSTRACT 

Real  earth  media  disperse  and  attenuate  propagating  waves.  This  anelastic 
behavior  can  be  well  described  by  a  viscoelastic  model.  We  have  developed  a  finite 
difference  simulator  to  model  wave  propagation  in  viscoelastic  media.  The  finite 
difference  method  was  chosen  in  favor  of  other  methods  for  several  reasons.  Finite 
difference  codes  are  more  portable  than,  for  instance,  pseudo- spectral  codes.  Moreover, 
finite  difference  schemes  provide  a  convenient  environment  to  define  complicated 
boundaries.  The  finite  difference  method  for  viscoelastic  wave  propagation  has  been 
thoroughly  investigated  with  convergence  tests,  dispersion  and  stability  analyses.  Two  2- 
D  examples  illustrate  the  difference  between  viscoelastic  and  elastic  wave  propagation. 
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INTRODUCTION 

To  a  high  degree,  elasticity  is  a  good  model  for  wave  propagation  through  the 
earth.  No  real  materials,  however,  are  perfectly  elastic,  but  rather  anelastic.  In  real 
media,  wave  energy  is  gradually  converted  into  heat.  Attenuation  of  propagating  wave 
forms  is  in  some  cases  quite  significant  and  could  be  a  source  of  erroneous  results  in 
forward  modeling,  inversion  and  imaging  if  neglected  ( Liu  et  al.,  1976).  The  quality 
factor,  Q,  characterizes  the  attenuation  of  waves  in  a  material.  This  is  a  measure  of  how 
many  wavelengths  a  wave  must  propagate  through  the  material  before  its  amplitude  drops 
by  1/e.  It  is  a  practical  measure  for  field  measurements  and  it  is  therefore  desirable  to 
develop  a  theoretical  model  using  quality  factors  for  pressure-  and  shear-waves. 

There  are  several  ways  of  modeling  anelastic  phenomena.  One  approach  is  to 
model  materials  through  their  mechanical  analogs  using  components  such  as  springs  and 
dashpots.  Such  a  model  is  called  a  viscoelastic  model  and  in  the  elastic  limit  this  reduces 
to  a  single  spring  with  a  spring  constant  corresponding  to  the  elastic  properties  of  the 
material.  The  so  called  standard  linear  solid  is  a  simple  viscoelastic  model  containing  of 
a  spring  in  parallel  with  spring  and  a  dashpot  in  series.  An  array  of  standard  linear  solids 
in  parallel  may  be  used  to  approximate  a  specific  viscoelastic  model  with  a  certain  Q 
versus  frequency  relation. 

Earth  materials  have  been  shown  to  have  a  nearly  constant  quality  factor,  Q,  over 
the  seismic  frequency  range  (review  in  Bourbie  et  al.,  1987).  A  viscoelastic  model 
consisting  of  a  series  of  standard  linear  solids  in  parallel  can  closely  approximate  a 
constant  Q  over  a  specified  frequency  range.  We  have  developed  a  constant  Q 
optimization  algorithm  based  on  work  by  Liu  et  al.  [1976].  As  a  rule  by  thumb  one 
standard  linear  solid  per  octave  in  the  specified  frequency  band  is  sufficient  to  achieve  a 
good  approximation.  Since  attenuation  and  dispersion  are  related  through  Kramers- 
Kronig's  relation,  the  constant  Q  model  also  yields  a  realistic  dispersion  relation. 
Futterman  [1962]  derived  a  theoretical  dispersion  relation  for  materials  with  a  constant 
Q.  Wuenschel  [1965]  experimentally  showed  Plexiglas  and  shale  to  have  nearly 
frequency-independent  quality  factor  over  several  octaves. 

We  based  our  finite  difference  approach  on  a  system  of  first  order  linear  partial 
differential  equations  derived  through  the  introduction  of  so  called  memory  variables. 
The  stress-strain  relation  contain  a  convolution  in  the  viscoelastic  case.  Viscoelastic 
media  may  therefore  be  regarded  as  remembering  earlier  deformation  states.  This  is 
reflected  through  the  memory  variables  which  are  added  to  the  equations  governing  the 
components  of  the  stress  tensor.  The  memory  variables  are  governed  by  separate 
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differential  equations.  Day  and  Minster  [1984]  also  derived  a  system  of  first  order  partial 
differential  equations  by  using  Pade  approximations  to  approximate  the  convolution  in 
the  stress-strain  relation. 

Control  of  numerical  dispersion  is  important  since  viscoelastic  media  are 
intrinsically  dispersive.  For  this  reason  earlier  viscoelastic  wave  propagation  simulators 
have  been  based  on  pseudo-spectral  methods  (e.g.  Carcione  et  al.,  1988).  Numerical 
dispersion  can  to  the  same  degree  also  be  controlled  using  finite  difference  methods  on 
denser  grids.  A  computationally  expensive  part  of  pseudo-spectral  schemes  is  the  FFT 
that  must  be  performed  each  time  step.  Moreover,  efficient  pseudo- spectral  codes  tend  to 
be  machine  dependent  due  to  FFT-routines,  etc.  This  is  not  the  case  for  finite  difference 
codes.  The  principal  arguments  for  choosing  a  finite  difference  method  to  simulate 
viscoelastic  wave  propagation  are  portability  of  efficient  code,  ability  to  model 
complicated  boundaries  accurately,  and  the  ability  to  control  the  simulation  accuracy. 

We  have  implemented  and  investigated  0(2,2),  0(2,4),  and  0(4,4)  accurate 
schemes.  These  have  been  thoroughly  investigated  through  convergence  tests,  dispersion 
and  stability  analyses.  The  discrete  dispersion  relation  for  the  system  was  found  to  be  the 
product  of  a  well  known  term  responsible  for  wave  propagation  and  another  term 
responsible  for  attenuation.  The  stability  criteria  for  the  viscoelastic  schemes  are 
approximately  the  same  as  for  the  analog  elastic  schemes.  The  Courant  number  must 
however  be  adjusted  to  the  highest  phase  velocity  (found  at  the  infinite  frequency)  in  the 
viscoelastic  medium.  We  have  performed  convergence  tests  where  we  compared  the 
analytic  solution  for  single  wavenumbers  to  the  numerical  solution  for  identical  initial 

conditions. 

In  this  paper  we  will  start  with  a  brief  review  of  viscoelastic  theory.  We  will  then 
present  our  algorithm  for  constant  Q  optimization.  We  will  then  pursue  the  theory  by 
deriving  the  equations  governing  viscoelastic  wave  propagation  and  also  formulate  these 
as  a  system  of  linear  first  order  partial  differential  equations.  Next  we  will  summarize  the 
implementation  and  the  tests  we  did  for  the  different  schemes.  After  this  we  will  present 
results  from  convergence  tests  and  stability  and  dispersion  analyses.  For  practical 
reasons  we  have  limited  this  analysis  to  1-D.  Our  results  may  be  directly  applied  to 
modeling  of  higher  dimensions  viscoelastic  wave  propagation.  We  describe  explicitly  the 
2-D  extension  of  our  method  and  exhibit  two  examples  where  disregarding  anelasticity 
seriously  limit  the  realism  of  the  problems. 
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THEORETICAL  MODEL  -  LINEAR  VISCOELASTICITY 

In  this  section  we  derive  a  theoretical  anelastic  model  based  on  viscoelastic 
theory,  a  phenomenological  way  to  describe  combined  elastic  and  viscous  behavior  of 
materials.  The  basic  hypothesis  is  that  a  current  value  of  the  stress  tensor  depends  upon 
the  history  of  the  strain  tensor.  The  simplest  viscoelastic  models  are  the  Maxwell  and  the 
Kelvin- Voigt  models,  which  essentially  represent  a  fluid  and  a  solid  respectively.  Both 
models  will  be  treated  more  thoroughly  below. 

Basic  definitions 

The  viscoelastic  hypothesis,  can  be  described  as. 


oo 


<rv(Q  =  ¥ij(ekl(t-s),ekl(t)) 

5=0 


(1) 


(, Christensen ,  1982)  where  y/  is  a  linear  tensor  valued  functional  which  transforms  each 
strain  history,  £<;/(r),-oo<f  <oo,  into  the  corresponding  stress  history,  <7- (r).  In  general 

all  field  variables  depend  both  on  time  and  position.  If  the  functional  is  linear  and  the 
strain  history  assumed  to  be  continuous,  Riesz’  representation  theorem  allows  us  to  write 
the  functional  as  a  Stieltjes  integral, 

tfy(0  =  je«(r-s)^w(s)-  (2) 

o 

The  fourth  order  tensor,  G,  collapses  into  two  independent  functions  for  an  isotropic  and 
homogeneous  material,  assuming  the  strain  history  to  be  a  symmetric  tensor  defined  as  in 
linear  elasticity.  By  including  a  step  discontinuity  for  the  strain  at  t=0,  constraining 
e  .(*)  =  0  for  t<0,  and  assuming  that  Gijkl  e  C1 ,  the  constitutive  relation  (2)  reduces  to, 

t  Ac  ( r) 

^(t)  =  Gm (0)£„ (r)  + 1 GijU ( t - 1)  — J dx  (3) 

0  WT 

In  1-D,  the  special  case  of  simple  shear  in  an  isotropic  homogeneous  material  (3)  reduces 
to, 
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CT(r)  =  G(0)e(0  +  \G{t-  t)e(r)dx 
0 


(4) 


where 


e(t)  = 


d£(t) 

dt 


G  is  called  the  stress  relaxation  function  and  is  the  viscoelastic  analog  to  the  Lame 
constant,  (I,  in  linear  elasticity.  Equation  (4)  may  also  be  written  as. 


<7=G*d£  =  £*dG  (5) 

where  *  denotes  Stieltjes  convolution.  For  sufficiently  well  behaving  functions  G  and  e, 
(5)  can  be  written  as, 


£7  =  G*  £=  £*G  (6) 

where  *  denotes  conventional  convolution. 

The  creep  compliance,  J,  is  the  inverse  of  G,  in  the  sense  that, 

£~a*j  =  J*a  0) 

By  Laplace  transforming  u_  and  (7)  and  eliminating  the  stress  and  the  strain,  we  obtain, 

sJ(s)sG(s)  =  1  (8) 

Viscoelastic  solids  and  fluids  may  be  distinguished  by  the  behavior  of  their  relaxation 
functions  at  infinite  time,  that  is  by  studying  G(°°).  If  G(°°)=0  the  material  is  defined  as  a 
fluid,  whereas  when  G^l^O  the  material  is  defined  as  a  solid. 

Differential  operator  form  of  stress-strain  relations  in  one  dimension 

If  we  assume  that  the  Laplace  transform  of  G  can  be  approximated  by  a  rational 
function,  this  is 
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sG(s)  = 


Q(s) 

P(s) 


then  the  stress-strain  relation  is  equivalent  to, 
P(D)a(t)  =  Q(D)e(t) 


(9) 


(10) 


where 


P(D)  =  ±pkDk 

k= 0 


do 


*= 0 


Here,  D  is  the  differential  operator  d/dt  and  N  is  an  arbitrary  natural  number.  This  is 
called  the  differential  operator  form  of  the  stress-strain  relations.  By  Laplace 
transforming  (10)  and  (11)  we  obtain. 


P(s)d(s)  --tPkts"°(k~r)(0)  =  <2U)e(s)  - -£<7*  i>r£(*"r)(0) 

S  k=\  r=  1  S  k= 1  r= 1 


(12) 


Through  identification  with  terms  in  (9)  and  (12)  we  obtain, 


'ZPr<jlr-‘\a)=Jiq,el'-l\0). 

r-k 


N 

L 

r=k 


(13) 


The  above  assumption  of  a  rational  approximation  can  be  justified  by  realizing  that  a 
realistic  relaxation  function  only  can  have  step  discontinuities,  which  results  in  a  1/s 
behavior  for  the  Laplace  transform  of  G  for  large  s.  Hence,  if  we  assume  the  order  of  Q 
to  be  less  than  or  equal  to  the  order  of  P+1  and  the  rational  sP(s)/Q(s)  to  have  no  multiple 
roots,  we  can  express  G  as  a  constant  plus  a  sum  of  N+l  (or  less)  inverse  first  order 
polynomials.  That  is, 


G(s)  =  K  + 


y  ak 

T?0bk  +5 


(14) 
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The  case  when  the  order  of  Q  is  greater  than  the  order  of  P+1  represents  infinite  stress  at 
finite  strain  or  in  other  words,  the  non-physical  case  of  rigid  bodies.  By  inverse  Laplace 
transforming  (14)  the  general  solution  for  G  in  the  time  domain  is  obtained  as, 

G(t)  =  K  +  fjcke~T'  (15) 

Jt=0 


A  mechanical  viscoelastic  model 

The  mechanical  models  used  in  viscoelasticity  are  usually  expressed  as 
combinations  of  elastic  and  viscous  elements.  The  elastic  model  is  the  spring  and  the 
viscous  model  is  the  dashpot  (illustrated  in  Figure  la).  The  constitutive  relations  for 
these  are, 

o  =  ke,  for  the  spring, 

and  (16) 

<j  =  rje,  for  the  dashpot. 

These  models  may  be  combined  either  in  series  or  in  parallel  to  form  different 
viscoelastic  elements.  A  spring  and  a  dashpot  connected  in  series  is  called  the  Maxwell 
model,  whereas  a  parallel  connection  is  called  the  Kelvin- Voigt  model.  The  Maxwell 
model  essentially  represents  a  fluid,  whereas  the  Kelvin- Voigt  model  represents  a  solid. 
Another  important  viscoelastic  configuration,  the  standard  linear  solid,  is  composed  of  a 
Maxwell  model  in  parallel  with  a  spring.  This  is  illustrated  in  Figure  lb.  A  similar 
configuration  is  the  so  called  Zener  model,  which  consists  of  a  Kelvin-Voigt  model  in 
series  with  a  spring.  Both  the  standard  linear  solid  and  the  Zener  model  are  viscoelastic 
solids. 

The  relaxation  function,  G,  for  the  standard  linear  solid  is, 

G(f)  =  kxe~“Ta  +  k2  (17> 


where 

and  the  creep  compliance,  J,  is, 
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K 

k 2(^1  "*■  kj ) 


e 


-t/Tc 


(18) 


where 


niK  +k2) 


Equation  (17)  may  also  be  written  on  the  form, 
G(t)  =  k2(  l  +  (^—\)e-lT°) 


L  standard  linear  solids  in  parallel  yield  the  relaxation  function,  G,  as, 

G(t)  =  ^(1-1(1-—)^) 

/=! 


(19) 


(20) 


where 

;=i 


ral  =  rll  /  ku 


*e'  =  V 


kl+K 

kxK 


Equations  (15)  and  (20)  are  equivalent  and  we  may  therefore  describe  an  essentially 
arbitrary  linear  viscoelastic  material  as  an  array  of  several  standard  linear  solids  in 
parallel.  This  is  also  equivalent  to  a  spring  in  parallel  with  an  array  of  several  Maxwell 
elements  in  parallel.  Expressing  the  causality  constraint  explicitly,  the  form  we  will  use 
for  a  general  viscoelastic  material  is. 


y{t)  =  MR 


f 

.  r 

\ 

\ 

1 

L 

-I 

1- 

IsL 

V 

/=i 

V 

/ 

e(t) 


(21) 
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where  Mp  is  the  relaxed  modulus  of  the  medium  ( Pipkin ,  1986)  and  6(t)  is  the  Heaviside 

A 

function  defined  as. 
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REALISTIC  ATTENUATION  AND  DISPERSION  MODELING 

Certain  experimental  results  have  shown  that  earth  materials  generally  have 
constant  Q,  both  for  pressure  and  shear  waves  over  a  limited  range  of  frequencies  (review 
in  Bourbie  et  al.,  1987).  Thus  it  is  highly  desirable  to  find  a  way  to  model  a  constant  Q, 
which  in  turn  will  yield  a  realistic  dispersion  relation  ( Futterman ,  1962).  In  this  section 
we  describe  an  algorithm  for  Q-optimization  for  an  array  of  standard  linear  solids  in 
parallel. 

Constant  Q  optimization 

The  quality  factor,  Q,  is  defined  as, 


Q(co)  = 


Re 

Mc(co) 

Im 

Mc(co) 

(22) 


where  Mc(cd)  is  the  complex  bulk  modulus,  defined  as  the  Fourier  transform  of  i g(t), 

which  for  the  case  of  L  standard  linear  solids  was  derived  in  (21). 

Our  objective  is  to  optimize  Q  such  that  it  is  approximately  constant  equal  to  <20 

over  a  limited  passband.  The  inverse  problem  is  better  formulated  if  we  optimize  1/Q 
instead  of  Q,  which  of  course  does  not  limit  us.  For  an  array  of  L  standard  linear  solids, 
this  is. 

Find 


Min  <p-  Min  J((?o  1  —Q  'f  dco 

Pass  Band 


where 


Id  ~  2 

l=\  \+(0  \ 


X-L^ 
1=1  1+' 


(23) 
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This  inverse  problem  is  ill  conditioned,  and  further  simplifications  are  necessary  to  obtain 

an  applicable  algorithm.  We  found  that  for  most  cases  of  interest  ( Bourbie  et  al.,  1987) 
one  can  make  the  approximation  t£[  =  tal .  This  approximation  and  (23)  yields. 


(=1  1  +  &rT, 


al 


(24) 


From  (23),  we  can  calculate  £T'  for  one  pair  of  relaxation  mechanisms.  This  is, 


Q- 1  =  (25) 

l  +  ft)2VT£ 

Note  that  when  rg/  =  xaV  (24)  can  be  regarded  as  the  sum  of  one-pair  relaxation 
mechanism  models  given  by  (25). 

For  one  pair  of  relaxation  mechanisms,  the  function  Q  has  a  unique  maximum 
at  a  certain  frequency  given  by, 


(O 


max 


1 

4^a 


(26) 


Since  rg/  =  we  use  the  notation  "to  fix  a  pair  of  relaxation  mechanisms  at  a  certain 
frequency"  when  setting, 


r 


Gl 


co 


(27) 


at  the  desired  frequency. 

Using  (23)  we  optimized  Q~x  with  respect  to  the  variables  (r^-r^)  and  ra(, 
using  the  conjugate  gradient  method.  In  Figure  2,  we  have  optimized  3  pairs  of 
relaxation  mechanisms  for  a  constant  Q  of  100,  between  100  and  500  Hz.  Note  that  the 

optimized  Q  is  very  close  to  the  desired  Q  in  the  frequency  interval.  It  should  be 
emphasized  that  the  approximation  r£l  =  is  better  the  larger  Q  is.  The  method  is 

slow,  but  as  we  shall  see  there  is  a  more  efficient  way  to  do  the  optimization.  In  Figure  3 
we  study  the  dependence  of  Q  on  perturbations  of  the  optimized  mechanisms.  The  effect 
of  increasing  all  relaxation  times  by  20%  is  minimal.  This  is  also  the  case  if  we  would 
increase  them  with  as  much  as  100%.  The  effect  is  much  more  dramatic  when  only  xal  is 
increased  by  20%,  while  keeping  xd  constant.  Thus,  a  small  perturbation  in  one  of  the 
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relaxation  times  in  a  pair  has  a  significant  impact  on  Q,  while  equal  perturbations  in  both 

relaxation  times  in  a  pair  have  negligible  effects.  From  the  experiments  we  draw  two 
conclusions.  First,  the  difference  between  r£l  and  r  ^  essentially  determines  the 

magnitude  of  Q.  Second,  when  optimizing  Q,  the  relaxation  times  will  be  evenly 
distributed  over  a  logarithmic  scale  ranging  from  the  lower  bound  of  the  pass  band  to  the 

upper  bound,  in  the  sense  discussed  above.  One  pair  of  relaxation  mechanisms  per  octave 
is  sufficient  to  yield  a  constant  Q  over  a  pass  band.  By  distributing  in  the  way 

described  above,  our  problem  is  reduced  to  an  over-determined  linear  system  of 
equations,  where  the  only  variables  are  (t^-  Tal).  This  was  solved  by  using  singular 

value  decomposition. 

In  Figure  4  we  have  optimized  3  relaxation  mechanisms  for  a  constant  Q  of  50, 
from  100  Flz  to  500  Hz.  Here  Q  is  overestimated  by  an  almost  constant  offset  (Q=52 
instead  of  Q=50).  The  error  is  larger  than  in  the  case  illustrated  above,  since  the 
approximations  are  worse  the  lower  Q  is.  This  is  due  to  the  difference  between  the 
relaxation  mechanisms  in  a  pair  which  increases  with  decreasing  Q.  If  we,  for  some 
reason,  want  to  have  a  very  accurate  approximation  of  a  constant  Q  of  50,  the 
optimization  may  be  done  at  Q=48  for  instance. 

We  now  have  an  algorithm  to  generate  a  realistic  attenuation  and  dispersion 
model  to  use  in  the  viscoelastic  code. 
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VISCOELASTIC  WAVE  PROPAGATION 

The  equations  describing  wave  propagation  in  viscoelastic  media  can  be  derived 
in  terms  of  the  creep  compliance  or  the  stress  relaxation  modulus.  Here  we  choose  to 
derive  it  from  the  latter.  For  simplicity,  we  will  derive  the  equations  for  the  one¬ 
dimensional  case  where  the  viscoelastic  equations  are  the  same  as  the  viscoacoustic.  A 
generalization  to  higher  dimensions  is  analogous  to  the  derivation  done  here.  The 
constitutive  equation  is, 


<j=  y/*S  (28) 

The  dot  denotes  a  time  derivative,  o  is  the  stress,  e  is  the  strain,  and  y  is  the  complex 
modulus  ( Pipkin ,  1986).  As  described  above  a  viscoelastic  medium  may  be  modeled  as 
an  array  of  standard  linear  solids  in  parallel.  Each  is  characterized  by  two  variables,  the 
stress  relaxation  time,  t0,  and  the  strain  relaxation  time,  r£.  The  relaxation  modulus  for 

an  array  of  L  parallel  standard  linear  solids  is  given  by  (21).  We  will  refer  to  a  particular 
set  of  relaxation  mechanisms  as  the  two  different  relaxation  times  corresponding  to  one 
of  the  steps. 

From  the  definitions  of  pressure  and  dilatation  ( Gurtin ,  1981)  we  know  that, 

(7  =  -p  (29) 

and 

e  =  vx  (30) 

Taking  the  time  derivative  of  (28)  and  using  (21)  leads  to, 

-p=  i j/*vx  (31) 


that  is, 
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The  convolution  terms  in  (32)  may  be  eliminated  by  introducing  so  called  memory 
variables,  which  will  be  denoted  rt  ( Carcione  et  al.,  1988).  Equation  (32)  thereby 

reduces  to, 
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From  (34)  we  see  that  r,  are  governed  by  convolutions  of  vx  with  exponential  functions, 
i.e.  the  kernels  of  r,  are  of  exponential  character.  A  set  of  first  order  linear  differential 
equations  may  be  obtained  instead  of  the  convolutions  as  follows.  First,  by  taking  the 
time  derivative  of  (34)  we  obtain, 
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which  gives  us, 
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From  (34)  we  see  that  (36)  reduces  to, 
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(37) 


We  have  now  derived  a  set  of  first  order  linear  differential  equations  for  the  memory 
variables.  Newton's  second  law  completes  the  full  description  of  wave  propagation  in  a 
viscoelastic  medium.  This  is, 


P*  =  -Px 


(38) 
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(33),  (37)  and,  (38)  is  the  system  of  first  order  linear  differential  equations  that  fully 
describe  one  dimensional  viscoelastic  wave  propagation  in  a  medium  with  L  sets  of 
standard  linear  solid  relaxation  mechanisms  that  we  will  use. 
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IMPLEMENTATION  AND  NUMERICAL  EXPERIMENTS 


The  system  of  L+2  first  order  linear  partial  differential  equations  we  wish  to  solve 
is  given  by  (33),  (37),  and  (38).  This  is, 
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We  made  several  tests  for  different  schemes  using  a  Ricker  wavelet  with  a  central 
frequency  at  35  Hz  propagating  in  the  positive  x  direction  in  a  2000  m  long  homogeneous 
interval  with  periodic  boundary  conditions.  The  initial  condition  is  plotted  in  Figure  5. 
The  relaxation  modulus  (corresponds  to  the  Lame  constant  in  the  elastic  case),  MR,  was  8 
GPa  and  the  density  2000  kg/m 3 .  This  leads  to  a  velocity  of  2000  m/s  in  the  0  Hz  limit. 
For  simplicity  we  used  one  set  of  relaxation  mechanisms  per  Q-model.  These  are  listed 
in  Table  1. 

Equation  (37)  proved  to  be  stiff  under  some  extreme  circumstances  but  otherwise 
quite  robust.  We  investigated  several  time  stepping  schemes  to  approximate  (37)  such  as, 
Adams-Bashforth  methods,  Gear's  methods,  and  the  Crank-Nicholson  scheme.  In  Figures 
6  and  7  we  display  two  cases  with  severe  numerical  dispersion  using  an  0(2,4)  scheme 
for  (33)  and  (38),  comparing  the  4:th  order  Gear's  method  to  Crank-Nicholson  as 
approximation  for  (37).  Even  in  these  extreme  cases  the  difference  is  minimal  and  hence 
we  decided  to  use  the  Crank-Nicholson  approximation  for  (37). 

We  decided  to  use  a  Leapfrog  scheme  to  approximate  the  time  derivatives  in  (33) 
and  (38).  This  conserves  energy  which  is  important  since  we  must  be  able  to  distinguish 
between  physical  and  numerical  attenuation  when  modeling  anelastic  media.  One  of  the 
principal  arguments  for  using  spectral  methods  instead  of  finite  differences  has  been  to 
avoid  confusing  physical  dispersion  with  numerical  dispersion.  If  all  spatial  derivatives 
in  our  system  of  differential  equations  is  approximated  by  second  order  accurate  central 
differences  the  numerical  dispersion  is  indeed  severe.  This  scheme  will  be  referred  to  as 
the  0(2,2)  scheme.  A  Courant  number  close  to  1  and  a  spatial  discretization  of  lm  (25 
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grid-points  per  wavelength  at  the  highest  frequency)  proved  to  be  optimal.  In  Figure  8 
we  compare  snapshots  at  2.25  s  for  the  0(2,2)  code  at  different  values  of  Q.  Q=10,000 
can  be  regarded  as  the  acoustic  limit  and  as  we  can  see  the  wavelet  has  suffered  from 
significant  numerical  dispersion.  As  we  shall  see  later,  this  scheme  also  does  not  model 
the  attenuation  correctly. 

Dablain  (1986)  and  Levander  (1988)  have  investigated  the  accuracy  of  high  order 
spatial  accurate  schemes  for  wave  propagation.  Increasing  the  spatial  accuracy  from 
second  to  fourth  order  improves  the  dispersion  characteristics  dramatically.  In  what  we 
will  refer  to  as  the  0(2,4)  scheme  all  second  order  accurate  approximations  of  the  spatial 
derivatives  have  been  replaced  by  fourth  order  accurate  approximations.  In  Figure  9  we 
compare  snapshots  at  2.25  s  at  different  values  of  Q.  There  is  no  visible  numerical 
dispersion  and  moreover,  as  we  shall  see,  the  attenuation  is  modeled  accurately.  In 
Figures  10  and  11  we  have  varied  the  Courant  number  and  the  spatial  discretization 
respectively.  A  Courant  number  of  0.2  and  a  spatial  discretization  of  2  m  (12  grid-points 
per  wavelength  at  the  highest  frequency)  proved  to  be  optimal.  This  scheme  has  several 
advantages.  First,  it  is  relatively  computationally  inexpensive  (compared  to  the  0(2,2) 
scheme  and  the  0(4,4)  scheme  which  we  will  describe  next).  It  also  models  wave 
propagation  in  viscoelastic  media  accurately  to  a  high  degree.  The  major  drawback  is  its 
poor  characteristics  at  high  Courant  numbers.  It  is  necessary  to  go  as  low  as  a  Courant 
number  of  0.2  in  order  to  have  sufficiently  low  numerical  dispersion.  On  the  other  hand, 
half  as  many  grid-points  per  wavelength  are  sufficient  compared  to  the  0(2,2)  scheme. 

The  0(2,4)  finite  difference  approximation  for  the  system  of  differential  equations 
is. 
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The  equivalent  differential  equations  (EDE)  are  derived  by  using  Taylor 
expansions  in  (39),  (40),  and  (41).  We  obtain. 
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Hence,  we  have  shown  that  the  schemes  are  second  order  accurate  in  time  and  fourth 
order  accurate  in  space. 

It  is  sometimes  possible  to  improve  the  accuracy  by  using  a  so  called  compact 
scheme  (Dab lain,  1986).  The  accuracy  of  finite  difference  approximations  such  as  (39) 
through  (41)  can  be  increased  from  second  to  fourth  order  in  time  by  subtracting  the 
highest  order  inaccurate  term  in  the  EDE  (equations  (42),  (43),  and  (44))  as  fourth  order 
spatial  finite  difference  approximations  using  the  exact  system  of  differential  equations 
(equations  (33),  (37),  and  (38)).  For  the  case  when  there  only  is  one  pair  of  relaxation 
mechanisms,  the  third  order  time  derivatives  expressed  as  spatial  derivatives  are. 
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We  used  (45)  and  (46)  in  (42)  and  (43)  to  derive  a  fourth  order  space  and  time  accurate 
compact  scheme  for  (33)  and  (38).  We  did  not  increase  the  accuracy  in  time  for  (37) 
since  this  did  not  improve  the  result  to  any  significant  extent.  We  will  refer  to  this 
scheme  as  the  0(4,4)  scheme. 

We  found  a  Courant  number  close  to  1  and  2m  grid-spacing  (12  grid-points  per 
wavelength  at  the  highest  frequency)  to  be  optimal.  This  is  illustrated  in  Figures  12  and 
13.  In  Figure  14  we  compare  snapshots  at  2.25  s  for  the  0(4,4)  scheme  at  different  values 
of  Q.  No  numerical  dispersion  is  visible.  The  accuracy  gained  compared  to  the  0(2,4) 
scheme  is  minimal  but  we  are  now  able  to  use  a  four  times  as  large  Courant  number. 

In  Figures  15  and  16  we  compare  snapshots  at  2.25  s  for  the  0(2,2),  the  0(2,4), 
and  the  0(4,4)  scheme  for  a  optimal  choice  of  simulation  parameters  at  Q=50  and 
Q=10,000.  The  solutions  from  the  0(4,4)  and  the  0(2,4)  schemes  are  close  with  minimal 
numerical  dispersion  while  the  0(2,2)  scheme  is  quite  poor. 

It  is  possible  to  obtain  analytic  viscoacoustic  solutions  for  some  simple 
geometries  by  replacing  the  velocity  in  the  Fourier  transform  of  the  acoustic  wave 
equation  by  the  complex  velocity  corresponding  to  the  specific  viscoacoustic  model 
(Bland,  1960).  This  allows  us  to  check  the  accuracy  of  our  finite  difference  schemes.  In 
Figures  17  and  18  we  compare  seismograms  using  two  different  values  of  Q  (Q=200  and 
Q= 10,000)  for  the  analytic  solution  and  the  0(4,4)  scheme.  Here  we  used  an  amplitude 
modulated  50  Hz  signal  as  our  initial  condition.  The  two  solutions  are  very  close  to  each 
other.  Since  we  already  have  showed  that  solutions  from  the  0(2,4)  and  the  0(4,4) 
scheme  are  close  we  have  now  showed  that  both  the  0(4,4)  scheme  and  the  0(2,4) 
scheme  model  viscoelasticity  very  accurately. 

In  Table  2  we  display  three  pairs  of  relaxation  mechanisms,  obtained  through  the 
optimization  algorithm  described  above,  which  yields  an  essentially  constant  Q  of  200 
between  20  and  80  Hz.  In  Figure  19  we  compare  a  snapshot  at  2.25  s  and  a  Q  of  200 
using  the  0(2,4)  scheme  for  one  pair  of  relaxation  mechanisms  (Table  1)  to  an  identical 
simulation,  but  for  a  case  when  we  use  three  pairs  of  relaxation  mechanisms  (Table  2). 
The  difference  is  small  and  hence  it  seems  like  the  optimization  of  a  constant  Q  is  not  as 
crucial  as  might  have  been  anticipated  in  simulating  viscoelastic  wave  propagation.  For 
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practical  purposes  it  might  even  be  sufficient  to  use  even  fewer  pairs  of  relaxation 
mechanisms  than  one  per  octave,  as  we  stated  as  good  rule  of  thumb.  The  explanation  for 
this  is  found  from  the  memory  variables.  In  Figure  20  we  have  plotted  the  memory 
variables  corresponding  to  the  different  relaxation  mechanisms  for  the  case  where  we 
used  three  pairs  to  model  a  constant  Q.  The  magnitude  of  the  curves  varies  dramatically. 
Hence,  the  main  contribution  to  the  viscoelastic  modeling  stems  from  the  memory 
variable  corresponding  to  the  first  set  of  relaxation  mechanisms  which  intuitively  seems 
to  be  reasonable.  In  Figure  21  we  compare  this  memory  variable  to  the  memory  variable 
in  the  simulation  where  we  only  used  one  pair  of  relaxation  mechanisms.  The  two  curves 
are  very  close  to  each  other. 

Finally,  in  Figure  22,  we  compare  two  snapshots  (t=1.25  s  and  t=2.25  s)  for  a 
simulation  using  the  three  relaxation  mechanisms  listed  in  Table  2  and  the  0(4,4)  code. 

We  found  that  the  0(4,4)  and  the  0(2,4)  scheme  model  viscoelasticity  sufficiently 
accurate  and  that  the  0(2,2)  scheme  hardly  is  of  any  use.  The  0(4,4)  scheme  is  only 
slightly  more  accurate  than  the  0(2,4)  scheme  for  moderate  Courant  numbers  and  is  also 
computationally  more  expensive.  This  is  reflected  through  the  number  of  calculations  per 
grid-point.  In  1-D,  the  number  of  calculations  per  grid-point  (counting  addition's, 
subtraction's,  multiplication's,  and  division's  equally),  is  about  75  for  the  0(4,4)  scheme, 
45  for  the  0(2,4)  scheme,  and  30  for  the  0(2,2)  scheme.  An  advantage  with  the  0(4,4) 
scheme  is  the  ability  of  running  simulations  at  high  Courant  numbers,  i.e.  fewer  time 
steps  are  needed,  which  makes  this  scheme  preferable  in  the  1-D  case.  To  achieve  high 
accuracy  over  a  wide  range  of  Courant  numbers  is  also  desirable  when  propagating  waves 
through  heterogeneous  media  where  velocities  may  vary  throughout  the  grid.  However, 
when  generalizing  to  the  2-D  viscoelastic  case,  the  0(4,4)  scheme  becomes  extremely 
computationally  expensive  so  that  the  advantages  with  the  0(2,4)  scheme  are  superior. 

The  pseudo  0(4,4)  scheme 

The  high  operation  count  and  the  asymmetry  of  the  dispersion  relation  for  the 
0(4,4)  scheme  led  us  to  derive  an  approximately  0(4,4)  scheme,  which  we  chose  to  call  a 
pseudo  0(4,4)  scheme.  The  viscoelastic  constitutive  relation  is  used  directly  instead  of 
the  equation  system  with  memory  variables.  Equation  (31)  and  (38)  are, 

P<=-Y*vx 

1  (47) 

v,= — Px 

l  P 
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To  create  a  compact  scheme  we  must  find  the  third  order  time  derivatives  of  pressure  and 
velocity  expressed  as  spatial  derivatives.  From  (47)  these  are  found  to  be 
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These  correction  terms  in  a  finite  difference  scheme  yields. 
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The  A  denotes  finite  difference  approximations  of  the  derivatives.  The  time  derivative  of 
the  relaxation  function  contains  one  term  multiplied  with  a  Dirac  delta  function  and  one 
term  multiplied  with  a  unit  step  function.  To  derive  the  pseudo  0(4,4)  scheme  the 
relaxation  function,  appearing  in  the  convolutions  with  the  third  order  spatial  derivatives, 
is  approximated  with  the  term  multiplied  with  the  Dirac  delta  function.  This 
approximation  yields. 
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Introduction  of  the  memory  variable  yields  the  following  system  of  equations,  for  one 
relaxation  mechanism, 
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Compared  to  the  0(4,4)  scheme  fewer  spatial  derivatives  have  to  be  calculated  in  each 
time  step.  The  accuracy  and  the  dispersive  properties  of  the  pseudo  0(4,4)  scheme  are 
close  to  those  of  the  0(4,4)  scheme,  as  is  demonstrated  in  the  Dispersion  and  Stability 
Analyses  section.  This  is  why  we  called  this  scheme  the  pseudo  0(4,4)  scheme. 
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CONVERGENCE  TESTS 

We  performed  a  series  of  convergence  tests  for  the  different  schemes  to  assure 
accuracy.  Convergence  was  shown  for  the  case  of  a  propagating  sinusoid,  since  an 
analytical  solution  for  this  problem  is  easy  to  obtain.  We  compared  this  to  the  numerical 
results  from  analog  finite  difference  simulations  and  calculated  the  difference  in  the  L2- 
sense. 

Analytic  solution  for  a  single  wavenumber 

An  analytic  solution  for  a  single  wavenumber,  k,  is  easy  to  find.  By  Fourier 
transforming  (33),  (37),  and  (38)  in  space  we  obtain  a  system  of  ordinary  differential 
equations  for  a  fixed  k. 


(52) 
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{-ik  Ip  0  0  , 

This  has  the  exponential  solution  which  is  obtained  through  the  eigenvalues  and 
corresponding  eigenvectors  of  A.  One  eigenvalue  is  real  and  corresponds  to  an 

exponentially  decaying  solution.  The  other  two  eigenvalues  are  complex  conjugates  and 
correspond  to  propagating  modes.  By  choosing  one  of  the  complex  eigenvalues,  the 
solution  for  a  sinusoid  propagating  in  one  direction  is  found.  The  characteristic  equation 
is  a  third  order  polynomial  equation  and  the  eigenvalues  and  eigenvectors  must  be 
determined  numerically.  We  chose  this  analytic  solution  for  convergence  tests  since 
these  can  b  determined  very  accurately. 


x-  Ax 


where 


x  = 


r 

\VJ 


A  — 


0 

0 
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Test  method 

The  convergence  tests  were  done  using  a  40  Hz  harmonic  wave  propagating  in  a 
medium  with  a  velocity  of  2000  m/s  and  Q=50.  The  initial  conditions  for  the  finite 
difference  schemes  were  calculated  using  the  analytic  solution.  Snapshots  after  1.00  s 
from  finite  difference  tests  were  compared  to  corresponding  analytic  solutions.  The  finite 
difference  schemes  were  tested  for  50,  25,  10,  5,  and  2  grid-points/wavelength.  The  ratio 
At/Ax  was  kept  constant  for  the  tests  using  the  0(2,2)  and  0(4,4)  schemes  so  that  the 
error  decreased  equally  in  time  and  space.  The  ratio  AtZ(Ax)2  was  kept  constant  when 
testing  the  0(2,4)  scheme  for  the  same  reason.  This  results  in  a  0(4,4)  behavior  for  the 
0(2,4)  scheme  during  the  convergence  tests.  The  0(2,4)  scheme  was  also  tested  for  a 
constant  Courant  number  of  0.4.  The  tests  were  performed  at  a  Courant  number  of  0.8 
for  the  0(2,2)  and  0(4,4)  schemes.  The  simulation  parameters  are  listed  in  Table  3. 

Results  of  the  convergence  tests 

The  difference  between  the  analytical  solution  and  the  results  using  the  0(2,2) 
scheme  is  displayed  in  Figure  23.  As  expected,  the  best  fit  line  for  the  logarithm  of  the 
L2-error  as  a  function  of  the  number  of  grid-points/wavelength  has  a  slope  close  to  -2  (- 
2.24).  The  surprisingly  close  fit  for  the  10  grid-points/wavelength  test  is  due  to  the 
periodicity  of  the  analytic  solution.  The  numerical  solution  is  actually  numerically 
dispersed  a  full  wavelength,  which  results  in  the  apparent  close  fit.  Snapshots  for 
different  numbers  of  grid-points/wavelength  at  1.00  s  are  plotted  in  Figures  24  and  25. 
The  result  for  the  0(2,4)  scheme  (non-constant  Courant  number)  is  displayed  in  Figure 
26.  The  best  fit  line  for  the  logarithmic  L2-error  has  a  slope  of  -5.13  which  is  close  to  the 
predicted  0(4,4)  behavior  (-4).  Snapshots  at  1.00  s  are  plotted  in  Figures  27  and  28.  The 
result  for  the  0(2,4)  scheme  using  a  constant  Courant  number  is  shown  in  Figure  29.  The 
best  fit  line  has  a  slope  of  -1.54.  Snapshots  at  1.00  s  are  plotted  in  Figures  30  and  31. 
The  result  of  the  0(4,4)  scheme  tests  is  presented  in  Figure  32.  The  best  fit  line  has  a 
slope  of  -4.23,  which  is  close  to  the  expected  (theoretically  it  should  be  -4).  The 
snapshots  at  1.00  s  are  plotted  in  Figures  33  and  34.  The  results  for  2  and  5  grid- 
points/wavelength  are  poor  for  all  the  schemes.  We  conclude  that  our  schemes  converge 
to  the  analytical  solution  as  predicted  by  theory. 
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DISPERSION  AND  STABILITY  ANALYSES 

Dispersion  and  stability  analyses  of  finite  difference  schemes  reveal  behavior  that 
otherwise  may  be  obscure  or  difficult  to  understand.  Viscoelastic  media  are  intrinsically 
dispersive.  The  phase  velocity  increases  with  frequency  as  shown  by  Futterman  [1962] 
and  Wuenschel  [1965].  This  is  evident  in  Figure  36.  In  addition  to  this  finite  difference 
methods  for  wave  propagation  introduce  numerical  dispersion  due  to  the  time  and  space 
discretization.  Hence,  we  are  studying  dispersive  media  with  dispersive  methods.  We 
must  therefore  assure  that  the  influence  of  numerical  dispersion  is  minimal. 

We  have  already  made  some  observations  of  dispersion  and  stability  issues  above 
under  Implementation  and  Numerical  Experiments.  A  more  thoroughly  study  is  however 
necessary  in  order  to  gain  deeper  understanding  and  to  confirm  our  initial  conclusions. 
For  simplicity,  in  the  following  we  will  consider  the  case  when  we  have  one  memory 
variable.  The  discrete  finite  difference  schemes  are  Fourier  (discrete)  transformed  to  the 
frequency-wavenumber  domain,  where  they  are  compared  to  the  analytical  dispersion  of 
the  viscoelastic  medium.  The  numerical  dispersion  originates  from  two  places  in  our 
schemes.  First,  unwanted  dispersion  might  be  introduced  in  the  equation  governing  the 
memory  variable,  r.  More  familiar  numerical  dispersion  is  introduced  through  the  part  of 
the  scheme  governing  the  wave  propagation.  The  dispersion  relation  for  the  viscoelastic 
schemes  must  be  solved  numerically,  most  convenient  through  an  iterative  method. 
Stability  criteria  were  investigated  using  von  Neumann  analysis.  The  numerical  scheme 
is  Fourier  transformed  in  space  and  Z  transformed  in  time.  The  roots  of  the  emerging 
polynomial  equation  yield  information  about  the  stability  of  the  scheme.  If  the  absolute 
value  of  all  roots  is  less  than  one  the  scheme  is  stable. 

Dispersion  analysis 

The  discrete  dispersion  relations  for  the  0(2,2)  and  0(2,4)  schemes  are  the 
product  of  a  term  responsible  for  wave-propagation  and  another  term  responsible  for 
viscoelastic  damping  and  dispersion,  originating  from  the  memory  variable  equation. 
The  relation  for  the  0(4,4)  scheme  does  not  have  the  same  form.  We  will  start  by 
investigating  the  equation  for  the  memory  variable,  r,  since  the  new  part  of  the  full 
dispersion  relation  is  coming  directly  from  that  equation. 

A  conventional  leapfrog  scheme  was  used  to  update  the  stress  and  velocity 
equations.  Since  a  dominating  part  of  the  numerical  dispersion  originate  from  the  wave 
propagating  part,  0(2,4)  and  0(4,4)  schemes  were  studied  in  addition  to  the  0(2,2) 
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scheme.  The  0(2,2)  and  0(2,4)  schemes  will  be  covered  in  the  same  paragraph,  since 
these  have  similar  dispersion  relations,  in  contrast  to  the  0(4,4)  scheme  which  is  a 
compact  scheme. 

Due  to  the  dependence  on  Q  and  the  complicated  analytical  dispersion  relations,  it 
is  not  possible  to  come  up  with  a  simple  rule  of  thumb  for  how  many  grid-points  per 
wavelength  that  is  needed  in  a  simulation  to  assure  sufficiently  low  effects  from 
numerical  dispersion. 


1.  The  memory  variable  equation. 

We  chose  to  use  an  unconditionally  stable  scheme  for  the  memory  variable,  since 
there  is  a  possibility  that  for  some  models,  the  stability  criteria  might  be  violated  for 
schemes  like  Euler  forward  and  Adams-Bashforths  method.  To  gain  higher  accuracy  we 
chose  the  Crank-Nicholson  above  the  Euler  backward  scheme.  The  damping  of  the  two 
schemes  are  compared  in  Figure  35,  where,  not  surprisingly,  the  Euler  backward  scheme 
does  not  yield  enough  damping  at  high  frequencies.  The  dispersion  relation  for  the 
Crank-Nicholson  scheme  is, 

4=l  +  (iM) - -  (53) 

®0  T(T  2  sin(twAt  /  2)  + — cos  (©A/  /  2) 


The  analytic  expression  for  the  dispersion  can  be  found  by  taking  the  limit  of  infinitely 
small  time  steps  At  ->  0 . 


Zr  1(0 

!  +  (—-!) - 


1(0- 


(54) 


where  cn()  is  the  angular  frequency  of  a  wave  propagating  in  an  elastic  medium. 


2.  The  0(2,2)  and  0(2,4)  schemes. 

The  dispersion  relations  for  the  0(2,2)  and  0(2,4)  schemes  may  be  written  as  a 
factor  for  the  numerical  dispersion  originating  from  an  elastic  scheme  multiplied  with  a 
factor  for  the  numerical  dispersion  from  the  equation  for  r.  The  dispersion  relations  for 
the  elastic  0(2,2)  and  0(2,4)  schemes  are, 
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sin2(£o0)  =  k2  sin2(/:0),  for  the  0(2,2)-scheme, 


(55) 


and 


sin2(<u0)  =  &2sin2(<:0) 


4-cos(^) 


-l2 


,  for  the  0(2,4)-scheme, 


(56) 


where  (O0  =  coAt,  kQ  =  kw Ax,  and  k  =  cAt/Ax  is  the  Courant  number.  Here  c  is  the  velocity 
of  the  medium  and  kw  is  the  wavenumber.  The  dispersion  relations  for  the  viscoelastic 
schemes  turns  out  to  be  as  simple  as  relation  (53)  multiplied  with  either  of  relation  (55)  or 
(56),  respectively.  Hence,  we  obtain, 


sin2(<y0)  =  k 2  sin2(£0)(l  +  F(co0,At,  xa,  te)) 


(57) 


sin2(ft)0)  =  k2  sin2  (>t0 ) 


4  -  cos  (Icq) 


(1  +  F((O0,At,  ta,Te)) 


(58) 


where  F  =  (<»2/ft)2)- 1  is  obtained  from  (53).  Equations  (57)  and  (58)  yield  that  the 
maximum  Courant  number  for  a  stable  scheme  is  lowered,  since  (1+F)  always  is  greater 
than  1.  The  stability  limit  (i.e.  the  maximum  Courant  number)  is  as  a  matter  of  fact 
decreased  by  approximately  a  factor  1/(1+1/Q),  since  F  is  almost  proportional  to  1/Q. 
The  highest  velocity  in  the  medium  is  found  at  infinite  frequency,  cmaj  r^p. 

Therefore,  to  determine  the  maximum  Courant  number  this  velocity  has  to  be  used.  The 
dispersion  curves  for  the  0(2,2)  and  0(2,4)  schemes  are  plotted  in  Figures  36  and  37. 
The  damping  curves  for  the  schemes  are  plotted  in  Figures  38  and  39.  The  viscoelastic 
dispersion  curves  demonstrate  the  same  characteristics  as  their  elastic  analogs.  To  model 
the  velocity  in  the  medium  with  sufficient  accuracy,  a  large  number  of  grid-points  per 
wavelength  is  needed  for  the  0(2,2)  scheme.  The  velocity  is  always  modeled  too  low. 
The  0(2,4)  scheme  overestimates  the  velocity  for  high  Courant  numbers  but  does  not,  on 
the  other  hand,  need  as  many  grid-points  per  wavelength  to  accurately  model  the  velocity. 
At  high  wavenumbers,  the  damping  is  underestimated  by  the  0(2,2)  scheme,  whereas  the 
0(2,4)  underestimates  the  damping  at  small  Courant  numbers  and  overestimates  the 
damping  at  large  Courant  numbers. 
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3.  The  0(4,4)  scheme. 


A  higher  order  scheme  as  the  0(4,4)  scheme  is  useful  since  it  models  wave 
propagation  accurately  for  a  wide  range  of  Courant  numbers.  This  is  an  important 
property  in  simulations  in  heterogeneous  media,  since  one  often  encounters  large  velocity 
variations  within  the  area  of  simulation.  The  compact  scheme  is  constructed  by  adding  a 
correction  term,  to  the  0(2,4)  scheme  consisting  of  spatial  derivatives,  as  was  described 
above.  Only  the  equations  for  stress  and  velocity  are  corrected  for  higher  time 
derivatives.  Hence,  the  scheme  used  for  the  memory  variable  equation  is  still  Crank- 
Nicholson.  SinCe  the  correction  is  not  done  through  higher  order  time  approximations  the 
dispersion  relation  is  no  longer  the  product  of  a  wave  propagation  term  and  a  relation 
from  the  memory  variable  equation.  The  dispersion  relation  weights  the  remainders  of 
the  0(2,4)  scheme  more  than  the  correction  terms.  For  low  Q,  this  results  in  a  behavior 
in  between  the  0(2,4)  and  0(4,4)  schemes.  The  dispersion  relation  for  the  0(4,4) 
compact  scheme  is, 


sin2(tu0)  = 


1 


1  +  i  t^2t;(/,Gt7  +  sin(*0))-^4 
6  sin(co0) 


f \  /4  ( 1  +  F)  +  \k2fff,(l  +  r,  )2  -  7*7(1  +  T,  )7 ]G/,/4  sin2  (fc0)  +  7  T]2t,/4(/,  rjG  +  sin  (k0 ) 


+  -  k2  T,  rj  sin2  (k0 )  sin  (co0 ) 
6 


(59) 


where 


F=  t,(1+  rjG) 


1 


cos(gj0  /  2) 


9  A  t 

i  sin(<y0  /  2)  -  - — cos (co0  /  2) 

2xc 


f.=efl+\k‘  <i+I,)/, 
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/,  =sin(£0) 


4-cos(fc0) 


/3  =  sin(*0)(cos(*0)-l) 


^ i  =  ~  1 


At 

T]  =  — 


At  the  elastic  limit  (  r,  — »  0)  the  expression  collapses  to  the  elastic  expression.  This  is. 


sin  2(a)0)  =  *J(/,+It2/,)2 


(60) 


which  is  equivalent  to, 


sin2(m0)  =  £2sin2(&0) 


^4-cos(£0)  +  &2(cos(fc0)-l7 
v  3  j 


(61) 


For  the  acoustic  scheme  the  stability  limit  is  k=1.0,  and,  as  mentioned  above,  slightly 
lower  for  the  viscoelastic  scheme.  Equation  (59),  yield  that  the  term  (1+F)  weights  f  ^ 

terms  more  than  f^  terms.  The  fj  terms  correspond  to  the  spatial  relation  for  the  0(2,4) 
scheme  while  the  f^  terms  correspond  to  the  correction  terms  in  the  0(4,4)  compact 

scheme.  For  low  Q  and  relatively  large  values  of  F,  this  leads  to  a  behavior  similar  to 
that  for  the  0(2,4)  scheme.  The  dispersion  curves  for  the  0(4,4)  scheme  are  plotted  in 
Figure  40.  The  fit  is  very  good  up  to  10  grid- points/wavelength.  The  difference  between 
the  theoretical  and  the  numerical  curve  is  less  than  1%.  The  damping  curves  for  the 
0(4,4)  scheme  are  displayed  in  Figure  41.  The  damping  is  overestimated  for  large  wave 
and  Courant  numbers. 

The  dispersion  relation  for  the  pseudo  0(4,4)  scheme  is,  as  for  the  0(2,2)  and 
0(2,4)  schemes,  almost  the  dispersion  relation  for  an  elastic  scheme  multiplied  with  a 
term  from  the  memory  variable  equation.  This  is, 

sin2(<y0)  =  k2{fx  +  ^2/3)2(  1  +  t,(1  +  r]G))  (62) 
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The  dispersion  curves  for  the  0(4,4)  and  pseudo  0(4,4)  schemes  are  compared  in  Figure 
42  and  the  damping  curves  are  compared  in  Figure  43.  The  dispersion  curves  are 
virtually  the  same  and  the  damping  behavior  of  the  pseudo  0(4,4)  scheme  is  actually 
slightly  better  than  for  the  0(4,4)  scheme.  Since  F  and  G  in  all  above  dispersion  relations 
depend  on  (Oo  in  a  sometimes  quite  complex  way,  co,  as  a  function  of  kw,  is  found  through 
iteration.  The  value  (D=ckw  serves  as  a  good  initial  value  in  this  procedure. 

Stability  analysis 

We  used  von  Neumann  analysis  to  investigate  stability  for  the  0(2,2),  0(2,4),  and 
the  0(4,4)  schemes.  This  method  consists  of  Fourier  transforming  the  spatial  difference 
approximations  of  the  finite  difference  scheme.  This  yields  a  difference  equation  for  each 
wavenumber.  If  the  absolute  values  of  the  solutions  to  the  difference  equation  are  less 
than  one  the  scheme  is  stable.  The  schemes  were  found  to  be  intrinsically  unstable.  The 
poles  of  the  modal  equation,  which  have  a  negative  real  part,  are  shifted  outside  the  unit 
circle,  whereas  the  poles  with  positive  real  part  are  shifted  inside.  This  should  be 
regarded  in  contrast  to  a  corresponding  elastic  scheme,  were  all  the  poles  are  on  the  unit 
circle.  In  Figure  44  we  have  plotted  the  poles  for  the  0(2,4)  scheme.  The  absolute  values 
of  the  poles  with  negative  real  part  are  plotted  in  Figure  45.  The  instability  remains,  even 
if  the  spatial  derivatives  are  computed  with  spectral  accuracy.  The  poles  for  this  pseudo- 
spectral  scheme  are  plotted  in  Figure  46.  Here,  the  poles  with  positive  real  part  are  found 
outside  the  unit  circle.  This  is  shown  in  Figure  47.  Increasing  the  Courant  number  or 
decreasing  Q  increases  the  absolute  values  of  the  poles  outside  the  unit  circle.  As  a 
consequence,  the  effect  of  the  instability  is  stronger  in  a  simulation  in  these  cases.  A 
simple  solution  to  overcome  this  problem  is  to  add  a  dissipative  term  in  one  or  both  of  the 
equations  for  the  pressure  and  velocity.  The  term  can  be  sufficiently  small  so  that  it  has 
no  effect  on  the  accuracy.  Equation  (63)  shows  an  example  how  to  add  the  dissipative 
term  in  the  equation  for  pressure. 

Pr=-— vx-r+* cpa  (63) 

T<r 

where  the  second  derivative  is  approximated  by, 


p^mAx)^ 


Pm+\  +  An-1  ~  2p, 
Ax2 


(64) 
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A  spatial  second  derivative  of  the  pressure  is  added  to  the  equation  to  yield  a  dissipative 
(heat-equation)  behavior.  The  parameter  k(x)  must  be  adjusted  in  accordance  with  the 
material  properties  of  every  grid-point,  k  has  to  be  large  enough  to  yield  a  stable  scheme 
and  still  not  affect  the  modeling  of  viscoelastic  waves.  The  second  derivative  must  be 
taken  in  the  (n-l):th  time  step  in  (39)  to  yield  the  desired  stabilizing  effect  Kreiss  and 
Oliger  [1973].  The  modal  equation  gives  a  few  hints  how  to  scale  k  to  achieve  a  stable 
scheme.  Presently  the  size  of  k  is  determined  by  a  formula  derived  from  the  modal 
equation  and  through  trial  and  error.  Figure  48  demonstrates  a  case  where  a  dissipative 
term  has  been  added  to  the  equation  for  the  stress  in  the  0(2,4)  scheme.  The  term 
controls  the  instability  as  is  illustrated  in  Figure  49,  where  the  absolute  values  of  the 
poles  are  less  than  or  equal  to  one.  The  effect  of  the  dissipative  term  on  dispersion  and 
damping  is  illustrated  for  the  0(2,4)  scheme  in  Figures  50  and  51,  where  these  are  plotted 
with  and  without  the  stability  controlling  dissipative  term.  At  high  wavenumbers  and  at 
the  peak  of  the  damping  curve,  the  damping  is  slightly  increased  by  the  dissipative  term. 
The  dispersion  curves  are  practically  identical. 

We  found  the  schemes  to  model  the  viscoelastic  equations  well.  The  propagation 
of  viscoelastic  waves  appears  to  be  the  most  crucial  issue,  since  the  Crank-Nicholson 
scheme  is  sufficiently  accurate  for  the  memory  variable  equation.  The  use  of  leapfrog 
time  stepping  schemes  seems  to  result  in  intrinsically  unstable  schemes.  Fortunately,  this 
problem  is  easy  to  overcome  by  adding  a  sufficiently  small  dissipative  term  which  does 
not  affect  the  solution.  We  are  now  ready  to  generalize  our  results  to  2-D. 
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2-D  VISCOELASTIC  MODELING 


The  theory  derived  for  the  1-D  case  above  is  easy  to  generalize  to  2-D  and  the 
results  from  the  1-D  experiments  are  applicable  also  in  the  2-D  case.  In  the  2-D 
viscoelastic  case  of  1  relaxation  mechanism  both  for  the  P-waves  and  the  SV-waves  the 
analog  to  (33),  (37),  and  (38)  is  a  system  of  9  first  order  linear  partial  differential 
equations.  This  is, 
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oi;  denotes  the  ij:th  component  of  the  symmetric  stress  tensor, 
v,  denotes  the  i:th  component  of  the  velocity. 

p,r,qvq2  are  the  memory-variables  (4  are  needed  for  the  simplest  2-D  case  with  one 

relaxation  mechanism  each  for  the  P-waves  and  the  SV-waves. 
xpr,  xp  are  the  strain  and  stress  relaxation  times  for  the  P-waves. 
x£,  rfa  are  the  strain  and  stress  relaxation  times  for  the  SV-waves. 

jU0  is  the  relaxation  modulus  corresponding  to  SV-waves  and  is  analog  to  the  Lame 
constant  /i  in  the  elastic  case. 

/r0  is  the  relaxation  modulus  corresponding  to  P-waves  and  is  analog  to  the  sum  of  the 
Lame  constant  A  and  2  times  the  Lame  constant  fi . 

Here  we  have  defined  the  relaxation  functions  so  that  these  correspond  to  P-  and 
SV-waves  respectively.  We  are  therefore  able  to  derive  the  relaxation  mechanisms 
directly  from  constraints  on  the  quality  factors  for  the  P-  and  SV-waves  separately, 
analogous  to  what  was  done  in  the  1-D  case  described  above.  This  is  why  a  relaxation 
modulus,  7T0 ,  not  corresponding  to  a  single  Lame  constant  is  introduced. 

In  the  simplest  1-D  case  when  only  one  pair  of  relaxation  mechanisms  was  used 
only  one  memory-variable  was  needed.  In  the  2-D  case  we  need  4  memory-variables  in 
the  simplest  case  corresponding  to  one  pair  of  relaxation  mechanisms  each  for  the  P-  and 
the  SV-waves.  This  stems  from  the  fact  that  different  convolution  terms  arises  in  the 
equations  governing  the  stress  components. 

Stability  and  accuracy  issues  are  treated  as  was  described  above.  The  scheme  has 
the  same  properties  as  the  1-D  0(2,4)  scheme.  Leap-frog  in  time  and  4:th  order  central 
differences  in  space  are  used  for  the  equations  governing  velocity  and  stress  components. 
Crank-Nicholson  in  time  and  4:th  order  central  differences  in  space  are  used  for  the 
equations  governing  the  memory  variables.  The  scheme  is  intrinsically  unstable  because 
it  yields  the  1-D  0(2,4)  scheme  when  applied  to  plane  waves  in  the  coordinate  directions, 
which  was  shown  to  be  unstable  above.  In  most  cases  the  instability  appears  to  be  very 
weak  and  the  scheme  can  be  used  directly  on  the  system  of  equations  in  (65). 

The  dispersion  relation  for  the  2-D  viscoelastic  scheme  has  two  solutions,  one 
corresponding  to  P-wave  velocity  and  one  corresponding  to  S-wave  velocity.  This  yields 
simple  modeling  of  fluid-solid  boundaries,  since  there  is  no  coupling  between  P-waves 
and  S-waves  in  a  homogeneous  medium.  Dissipative  terms  have  to  be  introduced,  as 
shown  above,  in  all  equations  for  stress  and  velocity  to  yield  a  conditionally  stable 
scheme  with  uncoupled  P-wave  and  S-wave  velocities,  when  this  is  necessary.  The 
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stabilizing  term  is  the  same  for  both  P-waves  and  S-waves,  though.  This  means  that  if 
the  S-waves  are  modeled  in  such  a  way  that  strong  dissipation  has  to  be  introduced  to 
keep  the  scheme  stable,  that  the  modeling  of  P-waves  might  be  affected. 

We  believe  that  our  viscoelastic  code  has  applications  in  a  diversity  of  fields.  To 
demonstrate  this  we  will  next  show  two  different  marine  environments,  one  continental 
shelf  and  one  deep  sea,  where  the  negligence  of  anelastic  effects  seriously  limits  realistic 
analyses.  In  these  cases  the  instability  does  not  pose  any  problems  and  no  stabilizing 
terms  were  therefore  added.  Close  to  the  mid-oceanic  ridges  the  seafloor  is  lineated 
parallel  to  the  ridges  ( Goff  and  Jordan,  1988).  The  geometries  of  the  environments  in  the 
two  examples  hence  have  symmetric  characters  along  one  axis,  and  a  2-D  study  is 
therefore  suitable  for  sufficiently  realistic  analyses. 

Numerical  example  1.  Realistic  deep-sea  seafloor  scattering. 

In  mono-static  seafloor  scattering  measurements  the  reverberation  of  acoustic 
energy  of  low  incidence  grazing  angle  on  the  seafloor  is  of  particular  interest.  Fluid 
saturated  and  porous  media,  such  as  sedimented  seafloor,  are  highly  anelastic.  This 
behavior  can  be  well  described  by  a  viscoelastic  model.  In  this  example  we  show  the 
importance  of  including  viscoelastic  effects  in  such  a  scattering  simulation.  The  seafloor 
and  the  parameters  used  are  typical  for  young  (10-50  m.y.  old)  seafloor  close  to  the 
northern  mid-Atlantic  ridge.  The  largest  scale  of  the  seafloor  obey  an  exponential 
decaying  solution  of  the  heat  equation  as  a  function  of  the  distance  to  the  mid  oceanic 
ridges.  When  modeling  the  seafloor,  larger  features  may  be  regarded  as  deterministic, 
while  a  stochastic  model  may  be  applied  to  smaller  features.  Goff  and  Jordan  [1988] 
present  a  new  way  of  analyzing  the  stochastic  component  of  the  seafloor,  assuming  it  is 
of  fractal  character.  Goff  et  al.  [1992]  describe  a  method  of  generating  realistic  stochastic 
seafloor  profiles  with  sediments  added.  We  used  a  profile  generated  by  them  from  data 
from  the  proximity  of  the  mid- Atlantic  ridge.  The  full  length  seafloor  profile  is  20  km 
long  and  is  discretized  down  to  1  m.  We  chose  a  1  km  segment  to  use  in  our  finite 
difference  simulations.  This  is  enclosed  in  Figure  52.  We  used  material  properties 
appropriate  for  seafloor  of  this  age  in  the  Atlantic  ocean  (e.g.  Hamilton  et  al.,  1982).  The 
sediments  are  characterized  by  a  low  Q  (30-50  both  for  pressure  and  shear  waves), 
pressure  wave  velocities  close  to  that  in  water,  low  shear  wave  velocities  (a  few  100  m/s), 
and  about  70%  higher  densities  than  that  of  water.  This  is  typical  for  fluid  saturated 
unconsolidated  sediments.  The  material  properties  are  listed  in  Table  4. 
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Our  initial  condition  was  a  tapered  amplitude  modulated  plane  wave  incident  on 
the  seafloor  with  10  degrees  grazing  angle.  The  wave  has  a  central  frequency  of  140  Hz 
and  a  bandwidth  of  approximately  50  Hz.  The  particular  wave-form  is  similar  to  wave¬ 
forms  used  in  many  sonar  systems.  Since  the  bandwidth  occupies  less  than  one  octave 
only  one  pair  of  relaxation  mechanisms  for  each  P-  and  SV-waves  is  needed.  The  plane 
wave  was  inserted  in  the  water  column  in  the  grid  as  an  initial  condition.  In  the  first  time 
step  the  desired  wave-front  was  chosen.  The  second  time  step  was  simply  calculated 
from  the  first  time  step  using  D'Alemberts  formula,  the  angle  of  incidence,  and  the 
velocity  in  the  water.  In  Figure  53a  we  show  the  P-wave  energy  snapshot  for  the  initial 
condition  (the  S-wave  energy  is  of  course  0  everywhere). 

We  chose  a  grid-spacing  of  1  m  in  both  spatial  directions.  The  grid  was  1 100  m 
times  315  m  with  an  additional  80  m  wide  frame  for  the  absorbing  boundary  added 
around  this  grid,  acting  as  a  sponge  filter.  In  this  frame  we  used  a  very  low  Q  (Q=6  both 
for  P-  and  S-waves)  and  tuned  the  Lame  constants  so  that  the  dispersion  relation  yielded  a 
constant  velocity  over  the  transition  from  the  grid  right  inside  the  frame  to  its  interior  for 
the  center  frequency.  Hence  we  were  able  to  minimize  reflections  from  impedance 
contrasts.  Below  we  will  present  a  suite  of  P-  and  SV-wave  energy  snapshots  .  We  have 
cut  away  the  absorbing  frame  in  all  of  these.  The  time  step  was  0.1  ms  which  leads  to  a 
maximum  Courant  number  of  0.26  which  is  well  below  the  stability  criterion  and  in 
accordance  with  the  accuracy  observations  for  the  0(2,4)  scheme.  It  should  be 
emphasized  that  if  we  sacrifice  accuracy  it  would  be  possible  to  run  simulations  at  a 
significantly  higher  Courant  number  if  a  radiating  boundary  condition  would  be  used 
instead  of  a  low  Q  frame  ( Higdon ,  1991),  since  the  regions  of  the  grid  with  the  highest 
Courant  number  are  within  the  absorbing  frame.  The  number  of  grid-points  per 
wavelength  varies  of  course  throughout  the  grid  both  with  P-  and  S-wave  velocities.  The 
lowest  number  of  grid-points  per  wavelength  occurs  for  SV-waves  in  the  sediments 
where  we  are  down  to  4,  while  it  is  19  for  the  P- waves  in  the  basement.  We  have  done 
convergence  tests  at  the  particular  material  properties  for  denser  sampled  grids  to  confirm 
that  the  field  is  not  under-sampled. 

One  has  to  keep  in  mind  that  a  viscoelastic  model  not  only  introduce  attenuative 
properties  to  the  materials  but  also  dispersive  effects  of  realistic  character  ( Futterman , 
1962).  In  this  example  we  can  therefore  expect  to  observe  two  different  effects  on  the 
scattered  field  when  comparing  elastic  and  a  viscoelastic  simulations.  Both  attenuation  of 
wave  energy  and  changes  in  impedance's  due  to  this  dispersion  will  affect  the  scattered 
field.  The  velocities  listed  in  Table  4  are  the  elastic  velocities  or  the  velocities  for  the  0 
Hz  limit  in  the  viscoelastic  case.  At  140  Hz  and  these  values  of  Q  the  velocity  increases 
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about  1%.  For  instance,  in  the  viscoelastic  case  at  140  Hz  the  velocity  for  the  P- waves  in 
the  sediment  layer  increases  from  1.54  to  1.55  km/s.  We  can  expect  the  effects  from 
these  changes  to  be  smaller  than  the  effects  from  the  attenuation,  although  these  might  be 
significant  if  we  are  close  to  transition  from  pre-  to  post-critical  reflection. 

In  Figures  53a  through  g  we  show  a  series  of  energy  snapshots.  Figure  53a  shows 
the  initial  condition.  We  present  snapshots  for  P-  and  SV-wave  energies  at  0.3  s  and  0.5  s 
in  Figures  53b  through  e.  At  0.3  s  there  is  a  strong  headwave  propagating  through  the 
basement  ahead  of  the  wave-front.  A  surfacewave  develops  in  the  interface  between 
sediments  and  the  basement.  At  0.5  s  this,  as  well  as  the  SV-wave  energy  in  the 
sediments,  has  been  attenuated  due  to  the  viscoelastic  properties  of  the  seafloor.  In 
Figures  53f  and  g  we  show  the  P-wave  energy  snapshot  at  0.5  s  for  an  elastic  sedimented 
seafloor  and  for  a  bare  seafloor.  The  surfacewave  is  much  stronger  in  the  elastic  case 
than  in  the  viscoelastic  case.  Hence,  the  viscoelastic  character  of  the  sediments  leads  to 
attenuation  of  the  surfacewave,  an  important  source  of  backscattered  energy  in  low 
frequency  elastic  scattering  simulations  (e.g.  Dougherty  and  Stephen,  1988).  The 
leftmost  crest  shadows  the  trough  in  the  bare  seafloor  experiment,  while  the  sediment 
layer  traps  energy  when  an  elastic  sediment  cover  is  present.  Therefore  the  surfacewave 
is  stronger  when  elastic  sediments  are  present  compared  to  the  bare  seafloor. 

We  have  done  several  different  experiments  with  numerous  seafloor  profiles 
which  all  show  that  anelastic  properties  play  an  important  role  in  this  kind  of  seafloor 
scattering  experiments.  It  is  beyond  the  scope  of  this  paper  to  go  into  the  details  of  these 
investigations  but  differences  up  to  15  dB  between  elastic  and  viscoelastic  simulations  are 
observed  in  the  back-scattered  direction. 

Numerical  example  2.  Wave  propagation  in  sediments  in  an  incised  valley  on  the 
continental  shelf. 

Common  features  on  the  continental  shelf  along  the  U.  S.  Gulf  coast  are  incised 
valleys  created  by  a  variety  of  fluvial  systems  (ancient  meandering  rivers  and  deltas,  etc.). 
During  retrogradation  such  systems  are  back-filled  with  clastic  sediments.  These 
sediment  fills  are  highly  anelastic  and  moreover  the  properties  changes  with  time  as  these 
are  compressed  by  younger  deposits.  We  chose  a  typical  cross-section  of  such  a  system 
for  a  numerical  simulation.  This  is  enclosed  in  Figure  54  and  the  material  properties  are 
listed  in  Table  5.  We  chose  not  to  model  the  water  surface  since  multiple  reflections  from 
this  would  complicate  a  simple  analysis  of  the  differences  between  a  viscoelastic  and  an 
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elastic  simulation.  Hence,  we  apply  the  same  type  of  absorbing  boundaries  all  around  the 
grid  as  was  used  in  the  first  example. 

As  a  source  we  used  an  80  Hz  Ricker  wavelet  modeling  an  air-gun  impulse.  The 
80  Hz  Ricker  wavelet  is  symmetric  with  respect  to  its  radius  and  was  inserted  as  a 
pressure  pulse  in  the  water  column,  constant  in  the  first  two  time  steps.  A  P-wave  energy 
snapshot  of  this  is  enclosed  in  Figure  55.  The  Ricker  wavelet  spans  over  3-4  octaves,  and 
hence  3  to  4  relaxation  mechanisms  should  be  applied  in  order  to  assure  a  constant  Q  as  a 
function  of  frequency  over  the  full  frequency  interval  to  be  in  accordance  with  the 
algorithm  discussed  above.  The  fewer  relaxation  mechanisms  that  are  used  the  less 
constant  Q  within  the  desired  interval.  It  should  however  be  emphasized  that  constant  Q 
not  is  a  universal  truth  for  all  earth  materials  but  deviations  from  this  rule  of  thumb  is  of 
course  possible.  The  purpose  of  this  example  is  to  qualitatively  demonstrate  the 
difference  between  an  elastic  and  viscoelastic  simulation  and  hence  we  chose  to  use  only 
one  set  of  relaxation  mechanisms.  This  leads  to  a  very  close  to  constant  Q  within  60  Hz 
to  120  Hz,  which  is  where  the  main  energy  content  of  this  Ricker  wavelet  is,  and  a  fairly 
good  approximation  within  the  rest  of  the  Ricker  wavelet  frequency  band.  This  serves 
our  purposes  more  than  well.  To  further  justify  our  decision  we  refer  to  the  1-D 
experiment  discussed  above,  which  is  illustrated  in  Figure  19,  where  we  compared 
snapshots  for  a  Ricker  wavelet  propagating  in  media  characterized  by  1  and  3  relaxation 
mechanisms,  respectively. 

The  model  is  190  m  in  the  vertical  direction  of  which  the  upper  70  m  are  water 
(see  Figure  54),  and  640  m  in  the  horizontal  direction.  We  have  also  in  this  example  cut 
out  the  absorbing  frame  in  all  energy  snapshots  in  Figure  55,  Figures  56a  through  h,  and 
Figures  57a  through  h.  1  m  grid-spacing  was  used  in  both  spatial  directions  and  the  time 
step  was  0.15  ms  which  leads  to  a  maximum  Courant  number  of  0.33  which  is  well  below 
the  stability  limit  for  our  scheme. 

In  Figures  56a  through  h  we  show  P-  and  SV-wave  energy  snapshots  of  a 
viscoelastic  simulation  using  the  material  properties  in  Table  5.  In  Figures  57  a  through  h 
we  show  snapshots  from  the  analog  experiment  in  the  elastic  limit  where  we  have  set  all 
Q  to  10,000  in  Table  5.  The  SV-wave  energy  snapshots  have  been  amplified  three  times 
more  compared  to  the  P-wave  energy  snapshots.  In  Figures  56a  and  b,  and  57a  and  b,  we 
show  the  P-  and  SV-wave  energy  snapshots  at  0.05  s  for  the  viscoelastic  and  the  elastic 
simulation  respectively.  In  the  snapshots  for  the  P-wave  energy,  reflections  from  the 
different  sediment  layers  are  clearly  visible.  There  is  a  strong  conversion  of  P-waves  to 
SV-waves  at  the  different  boundaries.  The  initiation  of  SV-waves  are  found  where  the  P- 
wave  intersect  the  boundaries  of  the  model.  This  gives  the  impression  that  the  SV-waves 
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are  following  the  boundaries  between  the  different  sediment  layers.  This  effect  is  not 
equally  prominent  at  the  water-solid  boundary.  There  is  hardly  any  difference  between 
the  elastic  and  viscoelastic  case.  The  energy  from  the  waves  is  also  well  absorbed  by  the 
boundaries  (no  reflections  are  visible).  Figures  56c  and  57c  are  snapshots  of  the  P-wave 
energy  after  0.1  s  of  propagation,  for  the  viscoelastic  respectively  the  elastic  model.  The 
P-wave  in  the  viscoelastic  medium  is  slightly  diminished  compared  to  the  P-wave  in  the 
elastic  medium.  The  conversion  of  P-waves  to  SV-waves  along  the  sediment  boundaries 
continues  as  can  be  seen  in  Figures  56d  and  57d,  which  are  snapshots  of  the  SV-wave 
energy  at  0.1  s  for  the  viscoelastic  respectively  the  elastic  medium.  The  SV-waves  in  the 
viscoelastic  medium  have  clearly  smaller  amplitudes  compared  to  the  SV-waves  in  the 
elastic  medium.  Details  are  also  more  difficult  to  observe  in  the  snapshot  from  the 
viscoelastic  simulation.  In  Figures  56e  and  f,  and  57e  and  f,  we  show  the  P-  and  SV- 
wave  energy  snapshots  at  0.2  s.  The  amplitude  of  the  P-headwave  is  significantly 
reduced  from  propagating  through  the  viscoelastic  sediments.  The  difference  between  the 
viscoelastic  and  the  elastic  simulations  is  even  more  striking  in  the  SV-wave  energy 
snapshots.  Moreover,  conversions  from  SV-to  P-waves  at  the  layer  boundaries  induce  P- 
wave  energy  in  the  water  column.  This  process  is  much  weaker  in  the  viscoelastic  case 
compared  to  the  elastic  case  (this  is  not  clearly  visible  in  the  snapshots  due  to  an 
amplitude  threshold  in  the  plotting  program).  In  Figures  56g  and  h,  and  57g  and  h,  we 
show  the  P-  and  SV-wave  energy  snapshots  at  0.3  s.  From  these  plots  it  is  evident  that 
wave  propagation  in  a  medium  like  this  is  of  viscoelastic  nature.  There  is  hardly  any 
energy  left  in  the  viscoelastic  simulation  compared  to  the  elastic  simulation. 
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CONCLUSIONS 

Viscoelasticity  provides  a  powerful  tool  for  modeling  real  earth  media.  The 
viscoelastic  relaxation  function  can  easily  be  tuned  to  resemble  the  attenuative  and 
dispersive  effects  that  real  earth  materials  have  on  propagating  waves.  We  conclude  that 
it  is  fully  possible  to  use  finite  difference  methods  to  simulate  wave  propagation  through 
dispersive  and  attenuative  viscoelastic  media.  It  is  easy  to  choose  simulation  accuracy 
and  the  code  is  highly  portable.  We  have  done  viscoelastic  wave  propagation  simulations 
using  our  schemes  and  have  achieved  high  performance  on  several  different  computers 
such  as,  the  Stardent  2000,  different  SUN  computers,  the  Cray  YMP,  and  the  hypercube 
i860. 

The  different  schemes  clearly  converge  towards  the  exact  solution  consistent  with 
predicted  accuracy.  There  are  two  different  instabilities  in  the  schemes.  First  there  is  the 
instability  from  the  finite  difference  approximation  of  the  elastic  wave  equation.  This  is 
possible  to  control  and  the  stability  limit  is  given  in  terms  of  the  Courant  number.  The 
other  instability  is  new  and  makes  the  schemes  intrinsically  unstable.  For  moderate 
simulation  parameters  and  Q  values,  this  tend  not  to  be  a  problem.  Generally  the  scheme 
is  more  unstable  the  lower  Q  is.  An  important  observation  is  that  low  values  of  Q  tend  to 
occur  where  the  velocity  is  low.  Since  the  highest  velocities  in  a  model  set  the  limit  for 
the  Courant  number,  the  Courant  number  in  low  velocity  zones  will  be  relatively  small 
and  therefore  also  less  unstable.  In  other  cases,  the  instability  can  be  controlled  by 
adding  a  sufficiently  small  dissipative  term  which  does  not  affect  the  accuracy  of  the 
scheme. 

The  0(2,4)  scheme  is  less  computationally  expensive  while  the  0(4,4)  scheme  has 
a  30%  higher  stability  limit  and  is  exact  over  a  significantly  larger  Courant  number 
interval,  which  is  important  in  wave  propagation  simulations  in  heterogeneous  media. 
Despite  the  advantages  of  the  0(4,4)  scheme,  the  computational  expense  in  2-D  makes  it 
impractical. 

Through  our  numerical  examples  we  have  shown  that  anelastic  effects  must  be 
considered  in  many  situations  and  that  our  scheme  accurately  models  wave  propagation 
through  viscoelastic  media.  Our  examples  also  prove  that  viscoelastic  simulations  using 
finite  difference  grids  of  practical  sizes  are  possible. 

The  code  has  primarily  been  developed  to  provide  a  useful  tool  for  scattering  and 
inversion  projects.  At  present  we  are  doing  2-D  viscoelastic  seafloor  scattering 
experiments  as  well  as  Q-inversion  studies. 
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Q 

xE(ms) 

%a(ms) 

50 

4.33 

4.16 

100 

4.287 

4.202 

200 

4.2654 

4.2230 

1,000 

4.2484 

4.2399 

10,000 

4.2446 

4.2437 

Table  1.  The  relaxation  mechanisms  used  in  the  experiments  to  simmulate  the  Q  in  the 
left  most  column. 
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Pair 

xE{ms) 

*  a(ms) 

1 

8.0186 

7.9578 

2 

1.1565 

1.1488 

3 

0.51773 

0.51705 

Table  2.  Optimized  relaxation  mechanisms  to  yield  a  constant  q  of  200  between  20  and 
80  Hz. 
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grid-points/wavelength  Ax  (m)  At  (ms) 

50  1  0.40 

25  2  0.80 

10  5  2.00 

5  10  4.00 

2  25  10.0 


Table  3a.  Simulation  parameters  for  the  0(2,2)  and  0(4,4)  scheme.  Q=50  @  40  Hz, 
velocity  in  medium  is  2.0  km/s  at  0  Hz. 


grid-points/wavelength  Ax  (m)  At  (ms) 

50  1  0.01 

25  2  0.04 

10  5  0.25 

5  10  1.00 

2  25  6.25 


Table  3b.  Simulation  parameters  for  the  0(2,4)  scheme.  Q=50  @  40  Hz,  velocity  in 
medium  is  2.0  km/s  at  0  Hz. 


grid-points/wavelentgth  Ax  (m)  At  (ms) 

50  1  0.20 

25  2  0.40 

10  5  1.00 

5  10  2.00 

2  25  5.00 


Table  3c.  Simulation  parameters  for  the  0(2,4)  scheme.  Q=50  @  40  Hz,  velocity  in 
medium  is  2.0  km/s  at  0  Hz.  (constant  Courant  number  =  0.4) 
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Water 

Sediments 

Basement 

vp  (km/s) 

1.54 

1.54 

2.6 

QP 

10,000 

50 

450 

vs  (km/s) 

0 

0.5 

1.45 

Qs 

0 

35 

225 

p  (g/cm3) 

1.05 

1.7 

2.8 

Table  4.  The  average  material  properties  used  in  the  seafloor  scattering  experiment.  vp 
is  the  velocity  of  pressure  waves  and  Qp  its  quality  factor,  vs  is  the  shear  wave  velocity 
and  Qs  its  quality  factor,  p  is  the  density. 
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Water 

Layer  1 

vp  (km/s) 

1.52 

1.60 

QP 

10,000 

40 

vs  (km/s) 

0 

0.4 

Qs 

0 

30 

p  (g/cm3) 

1.05 

1.3 

Layer  2 

Layer  3 

Layer  4 

1.75 

1.9 

2.2 

50 

50 

100 

0.8 

1.0 

1.2 

35 

45 

70 

1.5 

1.5 

2.0 

Table  5.  The  average  material  properties  used  for  the  ensized  valley  model.  vp  is  the 
velocity  of  pressure  waves  and  Qp  its  quality  factor,  vs  is  the  shear  wave  velocity  and 
Qs  its  quality  factor,  p  is  the  density.  Layer  1,  2,  3,  and  4  denotes  the  different  sediment 
layers. 
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FIGURE  CAPTIONS 

Figure  1.  Mechanical  models  of  viscoelastic  elements.  The  spring  and  the  dashpot  are 
illustrated  in  (a)  while  the  standard  linear  solid  is  illustrated  in  (b). 

Figure  2.  Three  optimized  pairs  of  relaxation  functions  for  a  constant  Q  of  100  between 
100  and  500  Hz.  Solid:  The  desired  function  (between  100  and  500  Hz).  Dashed:  The 
optimized  function  using  the  approximation  described  in  the  text.  Dotted:  Exact  Q  using 
the  optimized  relaxation  times. 

Figure  3.  Effects  on  Q  by  perturbations  in  the  relaxation  times.  The  case  illustrated  in 
Figure  1  was  used  as  the  unperturbed  state.  Solid:  The  desired  function  (between  100  and 
500  Hz).  Dashed:  Q  after  increasing  both  the  stress  relaxation  time  and  the  strain 
relaxation  time  by  20%.  Dotted:  Q  after  increasing  only  the  strain  relaxation  time  by 
20%. 


Figure  4.  Three  optimized  pairs  of  relaxation  functions  for  a  constant  Q  of  50  between 
100  and  500  Hz.  Solid:  The  desired  function  (between  100  and  500  Hz).  Dashed:  The 
optimized  function  using  the  approximation  described  in  the  text.  Dotted:  Exact  Q  using 
the  optimized  relaxation  times. 

Figure  5.  The  initial  condition  (the  stress)  in  the  numerical  experiments  is  a  Ricker 
wavelet  with  a  center  frequency  at  35  Hz,  traveling  in  a  2000  m  long  interval  with 
periodic  boundary  conditions  in  the  direction  of  increasing  x  with  a  velocity  of  2000  m/s. 
The  density  is  2000  kg  Inf* . 

Figure  6.  A  case  with  severe  numerical  dispersion  illustrating  robust  nature  of  the 
equation  governing  the  memory  variables.  The  same  initial  condition  and  material 
properties  as  in  Figure  5  with  a  Q  of  50  (relaxation  times  listed  in  Table  1)  were  used. 
The  snapshot  was  taken  at  2.25  s  using  a  Courant  number  of  0.6  and  a  step  size  of  2m  in 
the  0(2,4)  code.  Solid:  Crank-Nicholson  used  as  the  time  stepping  scheme  in  the 
equation  governing  the  memory  variables.  Dashed:  Fourth  order  Gear's  method  used  as 
the  time  stepping  scheme  in  the  equation  governing  the  memory  variables. 
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Figure  7.  A  case  with  severe  numerical  dispersion  illustrating  robust  nature  of  the 
equation  governing  the  memory  variables.  The  same  initial  condition  and  material 
properties  as  in  Figure  5  with  a  Q  of  10,000  (relaxation  times  listed  in  Table  1)  were 
used.  The  snapshot  was  taken  at  2.25  s  using  a  Courant  number  of  0.6  and  a  step  size  of 
2m  in  the  0(2,4)  code.  Solid:  Crank-Nicholson  used  as  the  time  stepping  scheme  in  the 
equation  governing  the  memory  variables.  Dashed:  Fourth  order  Gear's  method  used  as 
the  time  stepping  scheme  in  the  equation  governing  the  memory  variables. 

Figure  8.  Snapshots  at  2.25  s  illustrating  the  dependence  on  Q.  using  the  0(2,2)  code  at  a 
Courant  number  of  0.8  and  a  step  size  of  lm.  The  same  initial  condition  and  material 
properties  as  in  Figure  5  were  used.  The  relaxation  times  for  the  different  values  of  Q  are 
listed  in  Table  1.  Solid:  Q=50.  Dashed:  200.  Dotted:  10,000. 

Figure  9.  Snapshots  at  2.25  s  illustrating  the  dependence  on  Q  using  the  0(2,4)  code  at  a 
Courant  number  of  0.2  and  a  step  size  of  2m.  The  same  initial  condition  and  material 
properties  as  in  Figure  5  were  used.  The  relaxation  times  for  the  different  values  of  Q  are 
listed  in  Table  1.  Solid:  Q=50.  Dashed:  200.  Dotted:  10,000. 

Figure  10.  Snapshots  at  2.25  s  illustrating  the  dependence  on  the  Courant  number  using 
the  0(2,4)  code  at  a  step  size  of  2m.  The  same  initial  condition  and  material  properties  as 
in  Figure  5  and  a  Q  of  10,000  were  used.  The  relaxation  times  for  the  different  values  of 
Q  are  listed  in  Table  1.  Solid:  Courant  number  =  0.6.  Dashed:  Courant  number  =  0.4. 
Dotted:  Courant  number  =  0.2. 

Figure  11.  Snapshots  at  2.25  s  illustrating  the  dependence  on  the  spatial  step  size  using 
the  0(2,4)  code  at  a  Courant  number  of  0.2.  The  same  initial  condition  and  material 
properties  as  in  Figure  5  and  a  Q  of  10,000  were  used.  The  relaxation  times  for  the 
different  values  of  Q  are  listed  in  Table  1.  Solid:  Spatial  step  size  =  lm.  Dashed:  Spatial 
step  size  =  2m.  Dotted:  Spatial  step  size  =  3m. 

Figure  12.  Snapshot  at  2.25  s  illustrating  the  dependence  on  the  Courant  number  using 
the  0(4,4)  code  at  a  step  size  of  2m.  The  same  initial  condition  and  material  properties  as 
in  Figure  5  and  a  Q  of  10,000  were  used.  The  relaxation  times  for  the  different  values  of 
Q  are  listed  in  Table  1.  Solid:  Courant  number  =  0.8.  Dashed:  Courant  number  =  0.4. 
Dotted:  Courant  number  =  0.2. 
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Figure  13.  Snapshots  at  2.25  s  illustrating  the  dependence  on  the  spatial  step  size  using 
the  0(4,4)  code  at  a  Courant  number  of  0.8.  The  same  initial  condition  and  material 
properties  as  in  Figure  5  and  a  Q  of  10,000  were  used.  The  relaxation  times  for  the 
different  values  of  Q  are  listed  in  Table  1.  Solid:  Spatial  step  size  =  1m.  Dashed:  Spatial 
step  size  =  2m.  Dotted:  Spatial  step  size  =  3m. 

Figure  14.  Snapshots  at  2.25  s  illustrating  the  dependence  on  Q  using  the  0(4,4)  code  at 
a  Courant  number  of  0.8  and  a  step  size  of  lm.  The  same  initial  condition  and  material 
properties  as  in  Figure  5  were  used.  The  relaxation  times  for  the  different  values  of  Q  are 
listed  in  Table  1.  Solid:  Q=50.  Dashed:  200.  Dotted:  10,000. 

Figure  15.  Snapshots  at  2.25  s  comparing  the  different  schemes  to  each  other  using  Q  = 
50  (see  Table  1)  and  optimal  choices  of  simulation  parameters  for  the  different  schemes. 
The  same  initial  condition  and  material  properties  as  in  Figure  5  were  used.  Solid:  The 
0(4,4)  scheme  at  a  Courant  number  0.8  and  a  spatial  step  size  of  lm.  Dashed:  The 
0(2,4)  scheme  at  a  Courant  number  of  0.2  and  a  spatial  step  size  of  2m.  Dotted:  The 
0(2,2)  scheme  at  a  Courant  number  of  0.8  and  a  spatial  step  size  of  lm. 

Figure  16.  Snapshots  at  2.25  s  comparing  the  different  schemes  to  each  other  using  Q  = 
10,000  (see  Table  1)  and  optimal  choices  of  simulation  parameters  for  the  different 
schemes.  The  same  initial  condition  and  material  properties  as  in  Figure  5  were  used. 
Solid:  The  0(4,4)  scheme  at  a  Courant  number  0.8  and  a  spatial  step  size  of  lm.  Dashed: 
The  0(2,4)  scheme  at  a  Courant  number  of  0.2  and  a  spatial  step  size  of  2m.  Dotted:  The 
0(2,2)  scheme  at  a  Courant  number  of  0.8  and  a  spatial  step  size  of  lm. 

Figure  17.  Seismogram  at  x=500  between  2.0  and  2.5  s  comparing  the  0(4,4)  scheme  to 
an  analytical  solution  to  the  solution  from  the  0(4,4)  scheme.  Q  =  10,000  (see  Table  1) 
and  the  material  properties  are  the  same  as  in  Figure  5.  The  initial  condition  is  a  narrow 
band  amplitude  modulated  50  Hz  harmonic  wave.  Solid:  Analytical  solution  from  Bland 
[I960].  Dashed:  Solution  from  the  0(4,4)  scheme  using  a  Courant  number  of  0.8  and  a 
spatial  step  size  of  lm. 


51 


Viscoelastic  Finite  Difference  Modeling 


Figure  18.  Seismogram  at  x=500  between  2.0  and  2.5  s  comparing  the  0(4,4)  scheme  to 
an  analytical  solution  to  the  solution  from  the  0(4,4)  scheme.  Q  =  200  (see  Table  1)  and 
the  material  properties  are  the  same  as  in  Figure  5.  The  initial  condition  is  a  narrow  band 
amplitude  modulated  50  Hz  harmonic  wave.  Solid:  Analytical  solution  from  Bland 
[I960].  Dashed:  Solution  from  the  0(4,4)  scheme  using  a  Courant  number  of  0.8  and  a 
spatial  step  size  of  lm. 

Figure  19.  Snapshots  at  2.25  s  comparing  experiments  with  a  constant  Q  of  200  and 
three  relaxation  mechanisms  (listed  in  Table  2)  to  one  with  one  relaxation  mechanism 
(listed  in  Table  1)  using  the  0(2,4)  scheme  at  a  Courant  number  of  0.2  and  a  spatial  step 
size  of  2m.  The  same  initial  condition  and  material  properties  as  in  Figure  5  were  used. 
Solid:  Three  pairs  of  relaxation  mechanisms.  Dashed:  One  pair  of  relaxation 
mechanisms. 

Figure  20.  Snapshots  at  2.25  s  comparing  the  memory  variables  in  an  experiment  with  a 
constant  Q  of  200  and  three  relaxation  mechanisms  (listed  in  Table  2)  using  the  0(2,4) 
scheme  at  a  Courant  number  of  0.2  and  a  spatial  step  size  of  2m.  The  same  initial 
condition  and  material  properties  as  in  Figure  5  were  used.  Solid:  The  memory  variable 
corresponding  to  pair  1  of  the  relaxation  mechanisms  in  Table  2.  Dashed:  The  memory 
variable  corresponding  to  pair  2  of  the  relaxation  mechanisms  in  Table  2.  Dotted:  The 
memory  variable  corresponding  to  pair  3  of  the  relaxation  mechanisms  in  Table  2. 

Figure  21.  Snapshots  at  2.25  s  comparing  memory  variables  from  experiments  with  a 
constant  Q  of  200  and  three  relaxation  mechanisms  (listed  in  Table  2)  to  one  with  one 
relaxation  mechanism  (listed  in  Table  1)  using  the  0(2,4)  scheme  at  a  Courant  number  of 
0.2  and  a  spatial  step  size  of  2m.  The  same  initial  condition  and  material  properties  as  in 
Figure  5  were  used.  Solid:  The  memory  variable  corresponding  to  pair  1  of  the  relaxation 
mechanisms  in  Table  2  from  the  experiment  with  three  pairs  of  relaxation  mechanisms. 
Dashed:  The  memory  variable  from  the  experiment  using  one  pair  of  relaxation 
mechanisms  (listed  in  Table  1). 

Figure  22.  Snapshots  at  1.25  s  and  2.25  s  using  the  0(2,4)  code  at  a  Courant  number  of 
0.2,  a  step  size  of  2m.  The  same  initial  condition  and  material  properties  as  in  Figure  5 
were  used  with  a  constant  Q  of  200  using  3  pairs  of  relaxation  mechanisms  (listed  in 
Table  2).  Solid:  Snapshot  a  1.25  s.  Dashed:  Snapshot  at  2.25  s. 
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Figure  23.  Logarithmic  L2  error  for  the  0(2,2)  scheme  as  a  function  of  logarithm  of 
grid-points/wavelength.  The  best  fit  line  for  the  logarithmic  L2  error  decreases  as  -2.23, 
close  to  the  theoretical  value  -2.  Solid:  Best  fit  line.  Stars:  L2  error 

Figure  24.  Snapshot  of  propagating  sinusoid  after  1.00  s  for  the  0(2,2)  scheme  (velocity 
2.0  km/s  in  medium).  Solid:  Analytic  solution.  Dashed:  Numerical  result  5  grid- 
points/wavelength.  Dotted:  Numerical  result  2  grid-points/wavelength. 

Figure  25.  Snapshot  of  propagating  sinusoid  after  1.00  s  for  the  0(2,2)  scheme  (velocity 
2.0  km/s  in  medium).  Solid:  Analytic  solution.  Dashed:  Numerical  result  50  grid- 
points/wavelength.  Dotted:  Numerical  result  25  grid-points/wavelength.  Dash-dotted: 
Numerical  result  10  grid-points/wavelength. 

Figure  26.  Logarithmic  L2  error  for  the  0(2,4)  scheme  as  a  function  of  logarithm  of 
grid- points/wavelength.  The  best  fit  line  for  the  logarithmic  L2  error  decreases  as  -5.13, 
close  to  the  theoretical  value  -4.  Solid:  Best  fit  line.  Stars:  L2  error 

Figure  27.  Snapshot  of  propagating  sinusoid  after  1.00  s  for  the  0(2,4)  scheme  (velocity 
2.0  km/s  in  medium).  Solid:  Analytic  solution.  Dashed:  Numerical  result  5  grid- 
points/wavelength.  Dotted:  Numerical  result  2  grid-points/wavelength. 

Figure  28.  Snapshot  of  propagating  sinusoid  after  1.00  s  for  the  0(2,4)  scheme  (velocity 
2.0  km/s  in  medium).  Solid:  Analytic  solution.  Dashed:  Numerical  result  50  grid- 
points/wavelength.  Dotted:  Numerical  result  25  grid-points/wavelength.  Dash-dotted: 
Numerical  result  10  grid-points/wavelength. 

Figure  29.  Logarithmic  L2  error  for  the  0(2,4)  scheme  as  a  function  of  logarithm  of 
grid-points/wavelength  (constant  Courant  number=0.4).  The  best  fit  line  for  the 
logarithmic  L2  error  decreases  as  -1.54.  The  good  fit  at  low  number  of  grid- 
points/wavelength  depends  on  the  large  damping  of  the  sinusoid.  Solid:  Best  fit  line. 
Stars:  L2  error 

Figure  30.  Snapshot  of  propagating  sinusoid  after  1.00  s  for  the  0(2,4)  scheme  (constant 
Courant  number,  velocity  2.0  km/s  in  medium).  Solid:  Analytic  solution.  Dashed: 
Numerical  result  5  grid-points/wavelength.  Dotted:  Numerical  result  2  grid- 
points/wavelength. 
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Figure  31.  Snapshot  of  propagating  sinusoid  after  1.00  s  for  the  0(2,4)  scheme  (constant 
Courant  number,  velocity  2.0  km/s  in  medium).  Solid:  Analytic  solution.  Dashed: 
Numerical  result  50  grid-points/wavelength.  Dotted:  Numerical  result  25  grid- 
points/wavelength.  Dash-dotted:  Numerical  result  10  grid-points/wavelength. 

Figure  32.  Logarithmic  L2  error  for  the  0(4,4)  scheme  as  a  function  of  logarithm  of 
grid-points/wavelength.  The  best  fit  line  for  the  logarithmic  L2  error  decreases  as  -5.13, 
close  to  the  theoretical  value  -4.  Solid:  Best  fit  line.  Stars:  L2  error 

Figure  33.  Snapshot  of  propagating  sinusoid  after  1.00  s  for  the  0(4,4)  scheme  (velocity 
2.0  km/s  in  medium).  Solid:  Analytic  solution.  Dashed:  Numerical  result  5  grid- 
points/wavelength.  Dotted:  Numerical  result  2  grid-points/wavelength. 

Figure  34.  Snapshot  of  propagating  sinusoid  after  1.00  s  for  the  0(4,4)  scheme  (velocity 
2.0  km/s  in  medium).  Solid:  Analytic  solution.  Dashed:  Numerical  result  50  grid- 
points/wavelength.  Dotted:  Numerical  result  25  grid-points/wavelength.  Dash-dotted: 
Numerical  result  10  grid-points/wavelength. 

Figure  35.  Damping  for  Crank-Nicholson  and  Euler  backward  schemes  in  the  equation 
for  the  memory  variable.  The  curve  for  the  Crank-Nicholson  scheme  is  close  to  the 
theoretical  curve  at  all  frequencies.  The  Euler  backward  scheme  underestimates  the 
damping  at  high  wavenumbers.  Q=50  @  40  Hz,  At/xo=0.24.  This  ratio  between  the 
relaxation  time  and  the  time  step  is  very  large  and  will  seldom  be  encountered  due  to  the 
stability  criterion  and  dispersion  properties  for  a  full  scheme.  Solid:  Theoretical  curve. 
Dashed:  Crank-Nicholson.  Dotted:  Euler  backward. 

Figure  36.  Dispersion  for  the  0(2,2)  and  0(2,4)  schemes.  The  Courant  number 
increases  from  0.05,0.2  to  0.4  for  both  schemes.  The  0(2,4)  scheme  overestimates  the 
velocity  for  a  Courant  number  of  0.4.  The  0(2,4)  curve  approximates  the  theoretical 
curve  considerably  better  than  the  0(2,2)  curve.  Maximum  frequency  500  Hz.  Q=50  @ 
7.5  Hz.  Solid:  Theoretical  curve.  Dashed:  0(2,4).  Dotted:  0(2,2). 

Figure  37.  Dispersion  for  the  0(2,2)  and  0(2,4)  schemes.  The  Courant  number 
increases  from  0.6  to  0.7  for  both  schemes.  The  0(2,4)  scheme  overestimates  the 
velocity  for  these  Courant  numbers.  The  fit  of  the  0(2,2)  is  still  poor.  Maximum 
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frequency  500  Hz.  Q=50  @  37.5  Hz.  Solid:  Theoretical  curve.  Dashed:  0(2,4).  Dotted: 

0(2,2). 

Figure  38.  Damping  for  the  0(2,2)  scheme.  The  Courant  number  varies  from  0.4,  0.7  to 
0.9  .  The  damping  is  underestimated  for  all  Courant  numbers.  Maximum  frequency  500 
Hz.  Q=50@37.5Hz.  Solid:  Theoretical  curve.  Dashed:  0(2,2)  C=0.4  .  Dotted:  0(2,2) 
C=0.7  .  Dash-dotted:  0(2,2)  C=0.9  . 

Figure  39.  Damping  for  the  0(2,4)  scheme.  The  Courant  number  varies  from  0.05,  0.4 
to  0.7  .  The  damping  is  underestimated  for  Courant  number=0.05,  0.4  and  overestimated 
for  Courant  number=0.7  .  Maximum  frequency  500  Hz.  Q=50  @  37.5  Hz.  Solid: 
Theoretical  curve.  Dashed:  0(2,4)  C=0.05  .  Dotted:  0(2,4)  C=0.4  .  Dash-dotted:  0(2,4) 
C=0.7  . 


Figure  40.  Dispersion  for  the  0(4,4)  scheme.  The  Courant  number  varies  from  0.2,  0.4, 
0.6,  0.8  to  0.9  .  The  0(4,4)  aligns  very  well  to  the  theoretical  curve.  Maximum 
frequency  500  Hz.  Q=50  @  37.5  Hz.  Solid  (top):  Theoretical  curve.  Dashed  (bottom): 
0(4,4)  C=0.2  .  Dotted:  0(4,4)  C=0.4  .  Dash-dotted:  0(4,4)  C=0.6  .  Solid  (bottom): 
0(4,4)  C=0.8  .  Dashed  (top):  0(4,4)  C=0.9  . 

Figure  41.  Damping  for  the  0(4,4)  scheme.  The  damping  is  overestimated  for  large 
Courant  numbers  and  slightly  underestimated  for  Courant  number=0.2.  Solid  (bottom): 
Theoretical  curve.  Dashed  (bottom):  0(4,4)  C=0.2  .  Dotted:  0(4,4)  C=0.4  .  Dash- 
dotted:  0(4,4)  C=0.6  .  Solid  (top):  0(4,4)  C=0.8  .  Dashed  (top):  0(4,4)  C=0.9  . 

Figure  42.  Dispersion  for  the  0(4,4)  scheme  and  the  pseudo  0(4,4)  scheme  for  a  Courant 
number  =  0.8  .  The  dispersion  for  the  two  schemes  is  virtually  identical.  Solid: 
Theoretical  curve.  Dashed:  0(4,4)  scheme.  Dotted:  pseudo  0(4,4)  scheme. 

Figure  43.  Damping  for  the  0(4,4)  scheme  and  the  pseudo  0(4,4)  scheme  for  a  Courant 
number  =  0.8  .  The  pseudo  0(4,4)  scheme  overestimates  the  damping  slightly  at  the 
point  of  maximum  damping,  but  it  aligns  better  to  the  theoretical  curve,  for  high 
wavenumbers,  than  the  0(4,4)  scheme.  Solid:  Theoretical  curve.  Dashed:  0(4,4) 
scheme.  Dotted:  pseudo  0(4,4)  scheme. 
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Figure  44.  Poles  of  the  modal  equation  for  the  0(2,4)  scheme.  The  poles  with  negative 
real  part  lies  outside  the  unit  circle.  Poles  located  on  the  positive  real  axis  inside  the  unit 
circle  originates  from  the  memory  variable  equation.  Solid:  Unit  circle.  Stars:  Pole 
locations. 

Figure  45.  The  absolute  values  minus  one  for  the  poles  with  negative  real  part  for  the 
0(2,4)  scheme.  The  poles  corresponding  to  half  the  Nyquist  frequency  have  the  largest 
absolute  values.  Courant  number=0.4,  Q=50  @  37.5  Hz.  Solid:  Absolute  value  of  poles. 

Figure  46.  Poles  of  the  modal  equation  for  a  scheme  with  spectral  accuracy  for  spatial 
derivatives.  The  poles  with  positive  real  part  lies  outside  the  unit  circle.  Q=50  @  37.5 
Hz.  Solid:  Unit  circle.  Stars  and  crosses:  Pole  locations. 

Figure  47.  The  absolute  values  minus  one  for  the  poles  with  positive  real  part  for  the 
spectral  scheme.  Q=50  @  37.5  Hz.  Solid:  Absolute  values  of  poles. 

Figure  48.  Poles  of  the  modal  equation  for  the  0(2,4)  scheme.  Dissipation  introduced  in 
p  equation.  All  poles  are  on  or  inside  the  unit  circle.  Courant  number=0.4,  Q=50  @  37.5 
Hz.  Solid:  Unit  circle.  Stars:  pole  locations. 

Figure  49.  The  absolute  values  minus  one  for  the  poles  with  negative  real  part  for  the 
0(2,4)  scheme  after  dissipative  term  has  been  added  to  the  p  equation.  The  poles 
corresponding  to  half  the  Nyquist  frequency  have  the  smallest  absolute  values.  Courant 
number=0.4,  Q=50  @  37.5  Hz.  Solid:  Absolute  value  of  poles. 

Figure  50.  Dispersion  for  the  0(2,4)  scheme  with  and  without  dissipative  term.  The 
dispersion  curve  for  the  scheme  with  dissipative  term  completely  coincides  with  the 
dispersion  curve  for  the  scheme  without  a  dissipative  term.  Courant  number  =0.4,  Q=50 
@  37.5  Hz.  Solid:  Theoretical  curve.  Dashed:  Dispersion  curve  for  0(2,4)  with 
dissipative  term.  Dotted:  Dispersion  curve  for  0(2,4)  without  dissipative  term. 

Figure  51.  Damping  for  the  0(2,4)  scheme  with  and  without  dissipative  term.  The  curve 
for  the  0(2,4)  with  dissipative  overestimates  the  damping  for  high  wavenumbers  and  at 
the  peak  of  the  damping  curve.  Courant  number=0.4,  Q=50  @  37.5  Hz.  Solid: 
Theoretical  curve.  Dashed:  Damping  curve  for  0(2,4)  with  dissipative  term.  Dotted: 
Damping  curve  for  0(2,4)  without  dissipative  term. 


56 


Viscoelastic  Finite  Difference  Modeling 


Figure  52.  The  sedimented  seafloor  profile  used  in  the  seafloor  scattering  experiment 
(6.2  km-7.2  km  from  profile  v2602c,  Goff  et  al„  1992).  Solid:  Basement  profile. 
Dashed:  Sediment  profile. 

Figure  53.  Initial  condition  (a)  and  snapshots  of  the  scattered  energy  from  viscoelastic 
sedimented  (b-e),  elastic  sedimented  (f),  and  bare  seafloors  (g).  The  tapered  amplitude 
modulated  plane  wave  of  140  Hz  is  incident  at  10  degrees  grazing  angle.  The  material 
properties  are  listed  in  Table  4. 

The  absorbing  boundaries  are  cut  out  from  the  snapshots.  The  unit  of  the  scales  on  the 
axes  are  meters  and  are  relative  to  the  upper  left  comer  in  the  grid  when  the  absorbing 
boundaries  are  included. 

a.  Initial  condition,  P-wave  energy  (S-wave  energy  is  0). 

b.  P-wave  energy  for  a  viscoelastic  sedimented  seafloor  at  0.2  s. 

c.  S-wave  energy  for  a  viscoelastic  sedimented  seafloor  at  0.2  s. 

d.  P-wave  energy  for  a  viscoelastic  sedimented  seafloor  at  0.4  s. 

e.  S-wave  energy  for  a  viscoelastic  sedimented  seafloor  at  0.4  s. 

f.  P-wave  energy  for  an  elastic  sedimented  seafloor  at  0.4  s. 

g.  P-wave  energy  for  a  bare  seafloor  at  0.4  s. 

Figure  54.  The  model  used  in  the  incised  valley  experiment.  The  area  above  the  solid 
line  is  water.  The  area  between  the  solid  and  the  dashed  line  represents  Holocene 
transgressive  marine  deposits.  The  area  between  the  dashed  and  the  dotted  line  represents 
transgressive  marine  overlying  transgressive  estuarine  deposits.  The  area  between  the 
dotted  and  the  dash-dotted  line  represent  regressive  fluvial  deposits.  The  dotted  line 
indicates  the  bay  line,  the  boundary  between  the  regressive  and  transgressive  sequences. 
The  area  below  the  solid  line  and  above  the  dash-dotted  line  are  late  Pleistocene  to 
Holocene  deposits.  The  area  below  the  dash-dotted  line  represents  early  Wisconsinian 
regressive  deposits.  The  star  shows  the  position  of  the  source. 

Figure  55.  Initial  condition,  P-wave  energy  (S-wave  energy  is  0).  The  absorbing 
boundaries  are  cut  out  from  the  snapshots.  The  unit  of  the  scales  on  the  axes  are  meters 
and  are  relative  to  the  upper  left  comer  in  the  grid  when  the  absorbing  boundaries  are 
included. 
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Figure  56.  Snapshots  of  P-wave  energy  (a,  c,  e,  g)  and  SV-wave  energy  (b,  d,  f,  h)  for 
the  viscoelastic  incised  valley  model.  Material  properties  are  listed  in  Table  5.  The 
absorbing  boundaries  are  cut  out  from  the  snapshots.  The  unit  of  the  scales  on  the  axes 
are  meters  and  are  relative  to  the  upper  left  corner  in  the  grid  when  the  absorbing 
boundaries  are  included. 

a.  P-wave  energy  for  the  viscoelastic  model  at  0.05  s. 

b.  SV-wave  energy  for  the  viscoelastic  model  at  0.05  s. 

c.  P-wave  energy  for  the  viscoelastic  model  at  0.1  s. 

d.  SV-wave  energy  for  the  viscoelastic  model  at  0.1  s. 

e.  P-wave  energy  for  the  viscoelastic  model  at  0.2  s. 

f.  SV-wave  energy  for  the  viscoelastic  model  at  0.2  s. 

g.  P-wave  energy  for  the  viscoelastic  model  at  0.3  s. 

h.  SV-wave  energy  for  the  viscoelastic  model  at  0.3  s. 

Figure  57.  Snapshots  of  P-wave  energy  (a,  c,  e,  g)  and  SV-wave  energy  (b,  d,  f,  h)  for 
the  elastic  incised  valley  model.  Material  properties  are  listed  in  Table  5.  The  absorbing 
boundaries  are  cut  out  from  the  snapshots.  The  unit  of  the  scales  on  the  axes  are  meters 
and  are  relative  to  the  upper  left  comer  in  the  grid  when  the  absorbing  boundaries  are 
included. 

a.  P-wave  energy  for  the  elastic  model  at  0.05  s. 

b.  SV-wave  energy  for  the  elastic  model  at  0.05  s. 

c.  P-wave  energy  for  the  elastic  model  at  0. 1  s. 

d.  SV-wave  energy  for  the  elastic  model  at  0.1  s. 

e.  P-wave  energy  for  the  elastic  model  at  0.2  s. 

f.  SV-wave  energy  for  the  elastic  model  at  0.2  s. 

g.  P-wave  energy  for  the  elastic  model  at  0.3  s. 

h.  SV-wave  energy  for  the  elastic  model  at  0.3  s. 
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The  spring 


The  dashpot 

(a) 


Standard  linear  solid 


(b) 

Figure  1.  Illustration  of  the  spring  and  the  dashpot  (a)  and  the  standard  linear  solid  (b). 
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Figure  52. 
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